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This article reviews the current status of lattice-dynamical cal- 
culations in crystals, using density-functional perturbation theory, 
with emphasis on the plane-wave pseudo-potential method. Sev- 
eral specialized topics are treated, including the implementation 
for metals, the calculation of the response to macroscopic electric 
fields and their relevance to long wave-length vibrations in polar 
materials, the response to strain deformations, and higher-order 
responses. The success of this methodology is demonstrated with 
a number of applications existing in the literature. 
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I. INTRODUCTION 

The theory of lattice vibrations is one of the best estab- 
lished chapters of modern solid state physics and very few 
of the astonishing successes met by the latter could have 
been achieved without a strong foundation of the former. 
A wide variety of physical properties of solids depend 
on their lattice-dynamical behavior: infrared, Raman, 
and neutron-diffraction spectra; specific heats, thermal 
expansion, and heat conduction; phenomena related to 
the electron-phonon interaction such as the resistivity of 
metals, super-conductivity, and the temperature depen- 
dence of optical spectra, are just a few of them. As a 
matter of fact, their understanding in terms of phonons 
is considered to be one of the most convincing evidences 
that our current quantum picture of solids is correct. 

The basic theory of lattice vibrations dates back to 
the thirties and the treatise of Born and Huang (1954) is 
still considered today to be the reference textbook in this 
field. These early formulations were mainly concerned 
with establishing the general properties of the dynami- 
cal matrices — such as e.g. their symmetry and/or ana- 
lytical properties — without even considering their con- 
nections with the electronic properties which actually 
determine them. A systematic study of these connec- 
tions was not performed until the seventies (De Cicco and 
Johnson 1969, Pick, Cohen, and Martin 1970). The rela- 
tionship between the electronic and the lattice-dynamical 
properties of a system is important not only as a mat- 
ter of principle, but also because (and perhaps mainly 
because) it is only exploiting these relations that it is 
possible to calculate the lattice-dynamical properties of 
specific systems. 

The state of the art of theoretical condensed-matter 
physics and of computational materials science is such 
that it is nowadays possible to calculate specific proper- 
ties of specific (simple) materials using ab initio quantum 
mechanical techniques whose only input information is 
the chemical composition of the materials. In the specific 
case of lattice-dynamical properties, a large number of 
ab initio calculations based on the linear-response theory 
of lattice vibrations (De Cicco and Johnson 1969, Pick, 
Cohen, and Martin 1970) have been made possible 
over the past 10 years by the achievements of density- 
functional theory (Hohenberg and Kohn 1964, Kohn 
and Sham 1965), and by the development of density- 
functional perturbation theory (Zcin 1984, Baroni et 
al. 1987a) which is a method to apply the former within 
the general theoretical framework provided by the latter. 
Thanks to these theoretical and algorithmic advances, 
it is now possible to obtain accurate phonon dispersions 
on a fine grid of wave-vectors covering the entire Bril- 
louin zone, which compare directly with neutron diffrac- 
tion data, and from which several physical properties of 
the system (such as heat capacities, thermal expansion 



coefficients, temperature dependence of the band gap, 
and so on) can be calculated. 

The purpose of the present paper is to illustrate 
in some detail the theoretical framework of density- 
functional perturbation theory, including several techni- 
cal details which are useful for its implementation within 
the plane-wave pseudo-potential scheme. We will also 
provide a representative, though necessarily incomplete, 
choice of significant applications to the physics of insu- 
lators and metals, including their surfaces, alloys, and 
micro-structures. 



II. DENSITY-FUNCTIONAL PERTURBATION THEORY 

A. Lattice dynamics from electronic-structure theory 

The basic approximation which allows one to decouple 
the vibrational from the electronic degrees of freedom in 
a solid is the adiabatic approximation by Born and Op- 
penheimer (BO) (Born and Oppcnhcimcr 1927). Within 
this approximation, the lattice-dynamical properties of a 
system are determined by the eigenvalues £ and eigen- 
functions $ of the Schrodingcr equation: 



n 2 d 2 



E(H) U(R)=£$(R), (1) 



where R/ is the coordinate of the 7-th nucleus, Mi its 
mass, R = {R/} indicates the set of all the nuclear coor- 
dinates, and E(R) is the clamped-ion energy of the sys- 
tem, which is often referred to as the Born-Oppenheimer 
energy surface. In practice, E(R) is the ground-state en- 
ergy of a system of interacting electrons moving in the 
field of fixed nuclei, whose Hamiltonian — which acts onto 
the electronic variables and depends parametrically upon 
R — reads: 
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where Zj is the charge of the 7-th nucleus, — e is the elec- 
tron charge, and 7?at(R) is the electrostatic interaction 
between different nuclei: 
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The equilibrium geometry of the system is given by the 
condition that the forces acting on individual nuclei van- 
ish: 
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whereas the vibrational frequencies, to, are determined by 
the eigenvalues of the Hessian of the BO energy, scaled 
by the nuclear masses: 
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(5) 



The calculation of the equilibrium geometry and of 
the vibrational properties of a system thus amounts to 
computing the first and second derivatives of its BO en- 
ergy surface. The basic tool to accomplish this goal is 
the Hcllmann-Feynman (HF) theorem (Hcllmann 1937, 
Feynman 1939) which states that the first derivative of 
the eigenvalues of a Hamiltonian that depends on a pa- 
rameter A, H\, is given by the expectation value of the 
derivative of the Hamiltonian: 



= * 
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where is the eigenfunction of H\ corresponding to the 
E\ eigenvalue: H\^\ = E\^\. In the BO approxima- 
tion, nuclear coordinates act as parameters in the elec- 
tronic Hamiltonian, Eq. (2). The force acting on the 7-th 
nucleus in the electronic ground state is thus: 
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where >I>(r,R) is the electronic ground-state wave- 
function of the BO Hamiltonian. The BO Hamiltonian 
depends on R via the electron-ion interaction that cou- 
ples to the electronic degrees of freedom only through the 
electron charge density. The HF theorem states in this 
case that: 
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where Vr(i") is the electron- nucleus interaction, 
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and riR,(r) is the ground-state electron charge density cor- 
responding to the nuclear configuration R. The Hessian 
of the BO energy surface appearing in Eq. (5) is obtained 
by differentiating the HF forces with respect to nuclear 
coordinates: 
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Eq. (10) states that the calculation of the Hessian of 
the BO energy surfaces requires the calculation of the 
ground-state electron charge density, nn,(r), as well as 



of its linear response to a distortion of the nuclear ge- 
ometry, <9nR.(r)/<9R/. This fundamental result was first 
stated in the late sixties by De Cicco and Johnson (1969) 
and Pick, Cohen, and Martin (1970). The Hessian ma- 
trix is usually called the matrix of the inter- atomic force 
constants (IFC). 



B. Density Functional Theory 

According to the preceding discussion, the calculation 
of the derivatives of the BO energy surface with respect 
to the nuclear coordinates only requires the knowledge 
of the electronic charge-density distribution. This fact 
is a special case of a much more general property of the 
systems of interacting electrons, known as the Hohenberg 
and Kohn (1964) theorem. According to this theorem, no 
two different potentials acting on the electrons of a given 
system can give rise to a same ground-state electronic 
charge density. This property can be used in conjunction 
with the standard Raylcigh-Ritz variational principle of 
quantum mechanics to show that a universal functional, 1 
F[n(r)], of the electron charge density exists such that 
the functional: 



E[n(r)} = F[n(r)] + / n(r)V(r)dr 



(11) 



is minimized by the electron charge density of the ground 
state corresponding to the external potential V(r), under 
the constraint that the integral of n(r) equals the total 
number of electrons. Furthermore, the value of the mini- 
mum coincides with the ground-state energy. This theo- 
rem provides the foundation of what is currently known 
as density-functional theory (DFT) (Parr and Yang 1989, 
Dreizler and Gross 1990). It allows an enormous concep- 
tual simplification of the quantum-mechanical problem 
of the search of the ground-state properties of a system 
of interacting electrons, for it replaces the traditional de- 
scription based on wave-functions (which depend on 3N 
independent variables, N being the number of electrons) 
with a much more tractable description in terms of the 
charge density, which only depends on 3 variables. Two 
major problems hamper a straightforward application of 
this remarkably simple result: the fact that the form of 
the F functional is unknown, and that the conditions to 
be fulfilled for a function n(r) to be considered an accept- 
able ground-state charge distribution (and hence the do- 
main of the functional F) are poorly characterized. The 



1 By universal it is meant here that the functional is indepen- 
dent of the external potential acting on the electrons, though 
it obviously depends on the form of the electron-electron 
interaction. 
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second problem is hardly addressed at all, and one usu- 
ally contents oneself to impose the proper normalization 
of the charge density by the use of a Lagrange multi- 
plier. The first problem can be handled by mapping the 
system onto an auxiliary system of non interacting elec- 
trons (Kohn and Sham 1965) and by making appropriate 
approximations along the lines described in the next sub- 
section. 



1. The Kohn-Sham equations 

The Hohenberg and Kohn theorem states that all 
the physical properties of a system of interacting elec- 
trons are uniquely determined by its ground-state charge- 
density distribution. This property holds independently 
of the precise form of the electron-electron interaction. 
In particular, when the strength of the electron-electron 
interaction vanishes, F[n(r)] defines the ground-state ki- 
netic energy of a system of non-interacting electrons as 
a functional of its ground-state charge-density distribu- 
tion, Tn[n(r)]. This fact has been used by Kohn and 
Sham (1965) to map the problem of a system of interact- 
ing electrons onto an equivalent non-interacting problem. 
To this end, the unknown functional F[n(r)] is cast into 
the form: 



F[n(v)]=T [n(v) 



2 J |r-r'| 

+E xc [n(r)], (12) 



where the second term is the classical electrostatic self- 
interaction of the electron charge-density distribution, 
and the so-called exchange-correlation energy, E xc , is de- 
fined by Eq. (12). 2 Variation of the energy functional 
with respect to n(r) with the constraint that the num- 
ber of electrons is kept fixed leads formally to the same 
equation that would hold for a system of non-interacting 
electrons subject to an effective potential (also called self- 
consistent, SCF, potential) whose form is: 

V SCF (r) = V(r)+e 2 J J^Ltf + v xc (t), (13) 



where 



v xc {r) = 



SE XC 
Sn(r) 



(14) 



The exchange-correlation energy is the name we give to 
the part of the energy functional that we do not know how 
to calculate otherwise. For this reason, it has been named 
stupidity energy by Feynman (1972). Whether or not this is 
an useful concept depends on its magnitude with respect to 
the total functional and on the quality of the approximations 
one can find for it. 



is the functional derivative of the exchange-correlation 
energy, also called the exchange-correlation potential. 

The power of this trick is that, if one knew the effec- 
tive potential V SCF (r), the problem for non- interacting 
electrons could be trivially solved without knowing the 
form of the non- interacting kinetic-energy functional, To. 
To this end, one should simply solve the one-electron 
Schrodinger equation: 



2m dr 2 



+ V SCF (r) ip n (r) = e n ip n (r). 



(15) 



The ground-state charge-density distribution and non- 
interacting kinetic energy functional would then be given 
in terms of the auxiliary Kohn-Sham (KS) orbitals, 
V>n(r): 



N/2 

n(r) = 2^|V„(r)| 2 

71=1 

.2 N/2 



n=l 



(16) 
(17) 



where TV is the number of electrons, and the system is 
supposed to be non-magnetic, so that each one of the N/2 
lowest-lying orbital states accommodates 2 electrons of 
opposite spin. In periodic systems the index n running 
over occupied states can be thought as a double label, 
n = {v,k}, where v indicates the set of valence bands, 
and k is a wave-vector belonging to the first Brillouin 
zone (BZ). 

The ground-state energy given by Eqs. (11,12) can be 
equivalently expressed in terms of the KS eigenvalues: 

N / 2 2 r ( \ I i\ 

/ m „v~* e z / n(rm(r ) , , . 
E[n(r)} = 2 £ e n - - / ^'U rdr' 

n=l ' ' 

+E xc [n(r)] - J n(r)v xc (r)dr. (18) 

Eq. (15) has the form of a non-linear Schrodinger equa- 
tion whose potential depends on its own eigenfunctions 
through the electron charge-density distribution. Once 
an explicit form for the exchange-correlation energy is 
available, this equation can be solved in a self-consistent 
way using a variety of methods. 

2. Local-Density Approximation and beyond 

The Kohn-Sham scheme constitutes a practical way to 
implement DFT, provided an accurate and reasonably 
easy-to-use approximation is available for the exchange- 
correlation energy, E xc [n(r)]. In their original paper, 
Kohn and Sham (1965) proposed to assume that each 
small volume of the system (so small that the charge den- 
sity can be thought to be constant therein) contributes 
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the same exchange-correlation energy as an equal vol- 
ume of a homogeneous electron gas at the same density. 
With this assumption, the exchange-correlation energy 
functional and potential read: 



E xc [n(r)} 
v xc (n(r)) 



(n(r))n(r)dr, 

de xc (n) 



dn 



(19) 
(20) 



n=n(r) 



where e xc (n) is the exchange-correlation energy per parti- 
cle in an homogeneous electron gas at density n. This ap- 
proximation is known as the local density approximation 
(LDA). Approximate forms for e xc (n) have been known 
for a long time. Numerical results from nearly exact 
Monte-Carlo calculations for the homogeneous electron 
gas by Ceperley and Alder (1980) have been parameter- 
ized by Perdew and Zunger (1981) with a simple ana- 
lytical form. More accurate parameterizations have been 
recently proposed by Ortiz and Ballone (1994). All these 
different forms are very similar in the range of electron 
densities relevant to condensed-matter applications and 
yield very similar results. 

The LDA is exact in the limit of high density or of 
a slowly varying charge-density distribution (Kohn and 
Sham 1965). LDA has turned out to be much more suc- 
cessful than originally expected (see for instance Jones 
and Gunnarsson (1989)), in spite of its extreme simplic- 
ity. For weakly correlated materials, such as semiconduc- 
tors and simple metals, LDA accurately describes struc- 
tural and vibrational properties: the correct structure 
is usually found to have the lowest energy, while bond 
lengths, bulk moduli, and phonon frequencies are accu- 
rate within a few percent. 

LDA also has some well-known drawbacks. A large 
overestimate (~ 20%) of the crystal cohesive and molec- 
ular binding energies is possibly the worst failure of 
this approximation, together with its inability to prop- 
erly describe strongly correlated systems, such as e.g. 
transition-metal oxides. Much effort has been put in the 
search for better functionals than LDA (see for instance 
Perdew et al. (1999)). The use of gradient corrections 
(Becke 1988, Perdew et al. 1996) to LDA has become 
widespread in recent years. Gradient corrections are gen- 
erally found to improve the account of electron correla- 
tions in finite or semi- infinite systems, such as molecules 
or surfaces; less so in infinite solids. 

In general, DFT is a ground-state theory and KS eigen- 
values and eigenvector do not have a well defined physical 
meaning. Nevertheless, in the lack of better and equally 
general methods, KS eigenvalues are often used to es- 
timate excitation energies. The general features of the 
low-lying energy bands in solids obtained in this way are 
generally considered to be at least qualitatively correct, 
in spite of the fact that the LDA is known to substantially 
underestimate the optical gaps in insulators. 



C. Linear response 

In Sec. II. A we have seen that the electron-density 
linear response of a system determines the matrix of 
its IFC's, Eq. (10). Let us see now how this re- 
sponse can be obtained within DFT. The procedure de- 
scribed in the following is usually referred to as density- 
functional perturbation theory (DFPT) (Zein 1984, Ba- 
roni et al. 1987a, Gonze 1995b). 

In order to simplify the notation and make the argu- 
ment more general, we assume that the external potential 
acting on the electrons is a differentiable function of a set 
of parameters, A = {A^} (A^ = R/ in the case of lattice 
dynamics). According to the HF theorem, the first and 
second derivatives of the ground-state energy read: 



dE 
dX~ 

d 2 E 



dVx(r) 
d 2 Vx(r) 



n\(r)dr, 



n\(r)dr 



+ 



I 



On x (r) dV x (v) 
dXi dXj 



dr. 



(21) 



(22) 



The electron-density response, dn\(r)/dXi, appearing in 
Eq. (22) can be evaluated by linearizing Eqs. (16), (15), 
and (13) with respect to wave-function, density, and po- 
tential variations. Linearization of Eq. (16) leads to: 



N/2 



An(r) = 4Rc^CWA^(r), 



(23) 



n=l 



where the finite-difference operator A A is defined as: 



The super-script A has been omitted in Eq. (23), as well 
as in any subsequent formulas where such an omission 
does not give rise to ambiguities. Since the external po- 
tential (both unperturbed and perturbed) is real, each 
KS eigenfunction and its complex conjugate are degen- 
erate. As a consequence, the imaginary part of the sum 
appearing in Eq. (23) vanishes, so that the prescription 
to keep only the real part can be dropped. 

The variation of the KS orbitals, Aip n (r), is obtained 
by standard first-order perturbation theory (Messiah 
1962): 



where 



e„)|AV„> = —(AV SCF - Ae„)|V„>, (25) 



SCF ~ 2^dr^ + VsCF{r) 



(26) 
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is the unperturbed KS Hamiltonian, 

AV SCF (r)=AV(r)+e 2 J ^^dv' 

dv xc (n) 



+ 



dn 



An(r), (27) 



=n(r) 



is the first-order correction to the self-consistent poten- 
tial, and Ae„ = {i/j n \AV S cF\ipn) is the first order varia- 
tion of the KS eigenvalue, e n . 

In the atomic physics literature, an equation analo- 
gous to Eq. (25) is known as the Sternheimer equa- 
tion, after the work in which it was first used to cal- 
culate atomic polarizabilities (Sternheimer 1954). A self- 
consistent version of the Sternheimer equation was in- 
troduced by Mahan (1980) to calculate atomic polariz- 
abilities within DFT in the LDA. Similar methods are 
known in the quantum chemistry literature under the 
generic name of analytic evaluation of second- order en- 
ergy derivatives (Gerratt and Mills 1968, Amos 1987). In 
the specific context of the Hartree-Fock approximation, 
the resulting algorithm is called the coupled Hartree-Fock 
method (Gerratt and Mills 1968). 

Equations (23-27) form a set of self-consistent equa- 
tions for the perturbed system completely analogous to 
the KS equations in the unperturbed case — Eqs. (13), 
(15), and (16) — with the KS eigenvalue equation, Eq. 
(15), being replaced by the solution of a linear system, 
Eq. (25). In the present case, the self-consistency require- 
ment manifests itself in the dependence of the right-hand 
side upon the solution of the linear system. As AV scf (t) 
is a linear functional of Ara(r) which in turn depends lin- 
early on the Af/>'s, the whole self-consistent calculation 
can be cast in terms of a generalized linear problem. Note 
however that the right-hand side of Eq. (25) for Aip n de- 
pends through An on the solution of all the similar equa- 
tions holding for the Aip m (m ^ n). Hence, all the N 
equations, Eq. (25), are linearly coupled to each other, 
and the set of all the A^'s is the solution of a linear 
problem whose dimension is (ATM/2 x AM/2), M being 
the size of the basis set used to describe the ^>'s. The 
explicit form of this big linear equation can be worked 
out directly from Eqs. (23-27), or it can be equivalently 
derived from a variational principle, as explained in Sec. 
II. C. 3. Whether this big linear system is better solved 
directly by iterative methods or by the self-consistent so- 
lution of the smaller linear systems given by Eq. (25) is 
a matter of computational strategy. 

The first-order correction to a given eigen-function of 
the Schrodinger equation, given by Eq. (25), is often ex- 
pressed in terms of a sum over the spectrum of the un- 
perturbed Hamiltonian, 



AV-n(r) = V Vm(r) 



(ip m \AV SCF \ip n ) 



(28) 



running over all the states of the system, occupied and 
empty, with the exception of the state being considered, 
for which the energy denominator would vanish. Using 
Eq. (28), the electron charge-density response, Eq. (23), 
can be cast into the form: 



N/2 

n— 1 m^n 



(ip m \AV SOF \ip n ) 



(29) 



Eq. (29) shows that the contributions to the electron- 
density response coming from products of occupied states 
cancel each other, so that the m index can be thought 
as running onto conduction states only. This is equiv- 
alent to say that the electron-density distribution does 
not respond to a perturbation which only acts onto the 
occupied-state manifold (or, more generally, to the com- 
ponent of any perturbation which couples occupied states 
among each other). 

The explicit evaluation of Aip n (r) from Eq. (28) would 
require the knowledge of the full spectrum of the KS 
Hamiltonian and extensive summations over conduction 
bands. In Eq. (25), instead, only the knowledge of the 
occupied states of the system is needed to construct 
the right-hand side of the equation, and efficient iter- 
ative algorithms — such as conjugate gradient (Press et 
al. 1989, Stich et al. 1989, Payne et al. 1992) or mini- 
mal residual (Press et al. 1989, Saad and Schultz 1986) 
methods — can be used for the solution of the linear sys- 
tem. In this way the computational cost of the determi- 
nation of the density response to a single perturbation is 
of the same order as that needed for the calculation of 
the unperturbed ground-state density. 

The left-hand side of Eq. (25) is singular because the 
linear operator appearing therein has a null eigenvalue. 
However, we saw above that the response of the system 
to an external perturbation only depends on the compo- 
nent of the perturbation which couples the occupied-state 
manifold with the empty-state one. The projection onto 
the empty-state manifold of the first-order correction to 
occupied orbitals can be obtained from Eq. (25) by re- 
placing its right-hand side with — P c AV SCF \il)n) , where 
P c is the projector onto the empty-state manifold, and 
by adding to the linear operator on its left-hand side, 
H SCF — e n , a multiple of the projector onto the occupied- 
state manifold, P v , so as to make it non-singular: 



(H SCF + aP v - e n )\Aip n ) = -P c AV SCF \ip n ) 



(30) 



In practice, if the linear system is solved by the 
conjugate-gradient or any other iterative method and the 
trial solution is chosen orthogonal to the occupied-state 
manifold, orthogonality is maintained during iteration 
without having to care about the extra P v term on the 
left-hand side of Eq. (30). 

The above discussion applies to insulators where the 
gap is finite. In metals a finite density of states (DOS) 
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occurs at the Fermi energy, and a change in the orbital 
occupation number may occur upon the application of an 
infinitesimal perturbation. The modifications of DFPT 
needed to treat the linear response of metals have been 
discussed by de Gironcoli (1995) and will be presented in 
some detail in Sec. II. C. 4. 



1. Monochromatic perturbations 

One of the greatest advantages of DFPT — as compared 
to other non-perturbative methods to calculate the vi- 
brational properties of crystalline solids (such as, e.g. 
the frozen phonon or the molecular-dynamics spectral 
analysis methods) — is that within DFPT the responses 
to perturbations of different wave-lengths are decoupled. 
This feature allows one to calculate phonon frequencies 
at arbitrary wave-vectors q avoiding the use of super- 
cells and with a workload which is essentially indepen- 
dent of the phonon wave-length. To see this in some 
detail, we first rewrite Eq. (30) by explicitly indicating 
the wave- vector k and band index v of the unperturbed 
wave- function, iffi, and by projecting both sides of the 
equation over the manifold of states of wave- vector k + q. 
Translational invariance requires that the projector onto 
the k+q manifold, P k+q , commutes with H SCF and with 
the projectors onto the occupied- and empty-state man- 
ifolds, P v and P c . By indicating with P k+q P, = P k+q 
and P k+q P c = P k+q the projectors onto the occupied 
and empty states of wave-vector k + q, Eq. (30) can be 
rewritten as: 

(P SCF + a P k+q -6 k )|A^ k+q ) = 

-P k+q AV SCF |v£>, (31) 

where | Ai/> k+q ) = P k+q |Ai/> k ). By decomposing the per- 
turbing potential, AV SCF , into Fourier components, 



(32) 



where Aw q CF (r) is a lattice-periodic function, Eq. (31) 
can be cast into the form: 



k+q I 



k+qw k+q 



A« q c > k ), (33) 



where v' runs over the occupied states at k + q, w k 
and Au k+q are the periodic parts of the unperturbed 
wave-function and of the k + q Fourier component of its 
first-order correction, respectively, and the coordinate- 
representation kernel of the operator Pscf, ^scf (r, r') = 
(v\HscfW), is defined in terms of the kernel of the SCF 
hamiltonian, h^cF^, r ') = ( r \HscFW), by the relation: 



^ k + q (r, r') = c-^+^h° SCF (r, r / )S k +^ r ' . (34) 

Eq. (33) shows that the time-consuming step of the self- 
consistent process, Eq. (30), can be carried on working 
on lattice-periodic functions only, and the corresponding 
numerical workload is therefore independent of the wave- 
length of the perturbation. 

Let us now see how the other two steps of the self- 
consistent process, Eqs. (23) and (27), can be carried 
on in a similar way by treating each Fourier component 
of the perturbing potential and of the charge-density re- 
sponse independently. The Fourier components of any 
real function (such as An and An) with wave-vectors q 
and — q are complex conjugate of each other: An~ q (r) = 
(An q (r))*, and similarly for the potential. Because 
of time-reversal symmetry, a similar results applies to 
wave-functions: Au k+q (r) = (A\ k_q (r))*. Taking into 
account these relations, the Fourier component of the 
charge-density response at wave-vector q are obtained 
from Eq. (23): 



An q (r)=4V U k *(r)A W k+q (r). 



(35) 



Eq. (27) is a linear relation between the self-consistent 
variation of the potential and the variation of the electron 
charge-density distribution. The Fourier component of 
the self-consistent potential response reads: 

A, q CF (r) = A, q (r) + e 2 J ^le-^'W 

dV * c{n) An q (r). (36) 



+ 



dn 



n— n(r) 



The sampling of the BZ needed for the evaluation of 
(35) is analogous to that needed for the calculation of 
the unperturbed electron charge density, (16), and it re- 
quires in most cases an equal number of discrete k points. 
An exception to this rule occurs when calculating the re- 
sponse of insulators to macroscopic electric fields — as dis- 
cussed in Sec. II. C. 2 — and for the calculation of phonons 
in metals in the presence of Kohn anomalies — as dis- 
cussed in Sec. II. C. 4. 

In conclusion, Eqs. (33), (35), and (36) form a set of 
self-consistent relations for the charge-density and wave- 
function linear response to a perturbation of a wave- 
vector, q, which can be solved in terms of lattice-periodic 
functions only, and which is decoupled from all the other 
sets of similar equations holding for other Fourier com- 
ponents of a same perturbation. Thus, perturbations of 
different periodicity can be treated independently of each 
other with a numerical workload which is, for each per- 
turbation, of the same order as that needed for the un- 
perturbed system. 
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2. Homogeneous electric fields 

The electron density response to a homogeneous 
(macroscopic) electric field requires a special treatment. 
In fact, several electrostatic properties of an infinite 
solid are strictly speaking ill-defined in the long wave- 
length limit because the electrostatic potential describ- 
ing a homogeneous electric field E, (V^(r) = eE • r), is 

both unbound from below and incompatible with Born- 
von-Karman periodic boundary conditions. In the lin- 
ear regime, however, these pathologies can be efficiently 
treated in an elementary way. As in the harmonic ap- 
proximation the lattice-dynamical properties of polar in- 
sulators depend for long wave-lengths on the linear re- 
sponse to a homogeneous electric field (see Sec. II. D. 2), 
we will limit here our analysis to the linear regime only 
and postpone the discussion of non-linear electrostatic 
effects to Sec. II. F. 

From a mathematical point of view, the main difficulty 
in treating macroscopic electric fields stems from the fact 
that the position operator, r, is ill-defined in a periodic 
system and so are its matrix elements between wave- 
functions satisfying Born-von-Karman boundary condi- 
tions. The wave- function response to a given pertur- 
bation, Eq. (28), however, only depends on the off- 
diagonal matrix elements of the perturbing potential be- 
tween eigen-functions of the unperturbed Hamiltonian. 
Such matrix elements are indeed well defined even for a 
macroscopic electric field, as it can be seen by re- writing 
them in terms of the commutator between r and the un- 
perturbed Hamiltonian, which is a lattice-periodic oper- 
ator: 

(1> m \r\1, n ) = MHso " YUn \ Vm^n. (37) 

If the self-consistent potential acting on the electrons is 
local, the above commutator is simply proportional to 
the momentum operator: 



n 2 d 

? ' r = 7T- 

m or 



(38) 



Otherwise, the commutator will contain an explicit con- 
tribution from the non-local part of the potential (Baroni 
and Resta 1986b, Baroni et al. 1987a, Hybertsen and 
Louie 1987). 

When calculating the response of a crystal to an 
applied electric field, En, one must consider that the 
screened field acting on the electrons is: 



4ttP, 



(39) 



where P is the electronic polarization linearly induced by 
the screened (i.e. self-consistent) field, E: 



P = ~ / rA E n(r)dr. 

V Jv 



(40) 



In addition to the macroscopic screening, expressed by 
Eq. (39), the density response to an external macroscopic 
electric field, E , also involves microscopic Fourier com- 
ponents, 



N/2 

A E n(r)=4^C(r)A E ^(r), 



(41) 



71=1 



responsible for the so called local fields, that must be 
taken into account in the self-consistent procedure. 

Eq. (40) is of course well defined for any finite sys- 
tem. The electric polarization of a macroscopic piece of 
matter, however, is ill defined in that it depends on the 
details of the charge distribution at the surface of the 
sample. Nevertheless, the polarization linearly induced 
by a given perturbation is well defined, and Eq. (40) 
can in fact be recast into a boundary-insensitive form 
(Littlewood 1980). To see this, we use Eq. (23) and we 
obtain from Eq. (40): 



4e 
V 



N/2 



P„ = ^5Z<^„|r a |A E V» n ) 

n=l 

N/2 oo 



<W[H " r ' rJII ".«, m |A E *,>, 



n=l m=AT/2+l 



(e n Cm) 



(42) 



where the subscript a indicates Cartesian components. 
Let us introduce the wave- function V>"(r) defined as: 

CW = E ^(r) <M^M , (43 ) 



ip satisfies an equation of the same kind as Eq. (30) with 
the perturbing potential on its right-hand side replaced 

by [H SCF ,r a ]: 



(H SCF - e n )\ip") = P c [H SCF ,r a ]\ip n ). 



(44) 



The induced polarization, Eq. (42), can be recast into 
the form: 



4e 
V 



N/2 



]T(C|A E Vn>, 



(45) 



n=l 



where the anti-hermitian character of the commutator 
has been used. 

The first-order correction to a crystal wave-function 
due to a perturbing homogeneous electric field, E, is given 
by the response to the full, macroscopic and microscopic, 
perturbation: 



(H SCF -e„)|A E ^„) = -e^E a |C) 



-P c AV'f\il> n ), (46) 
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where: 



+ 



dv xc (n) 



dn 



A E n(r). 



(47) 



n— n(r) 



The terms of the sum appearing in (45) implicitly be- 
have as the inverse of the difference between a conduc- 
tion and a valence energy eigenvalue to the third power, 
~ (e c — e v )~ 3 . One of the energy denominators results 
from the standard first-order perturbation calculation of 
the perturbed wave function, A^ip ni Eq. (28). A sec- 
ond energy denominator comes from the term ip in the 
right-hand side of Eq. (46), which requires by itself the 
solution of a linear equation analogous to that of first- 
order perturbation theory, Eq. (44). The third energy 
denominator comes from the ip, which appears explicitly 
in the bracket of Eq. (45). This dependence of the terms 
of the sum upon the direct gaps at different k points of 
the BZ may require a rather fine sampling of the BZ when 
the fundamental gap is small. In these cases, the number 
of k points needed to compute the dielectric constant is 
substantially larger than that needed by a standard un- 
perturbed calculation (Baroni and Resta 1986a, Baroni 
and Resta 1986b, de Gironcoli et al. 1989). 

The self-consistent cycle defined by Eqs. (39-47) can 
be performed starting from a given external macroscopic 
electric field, E , and updating at each iteration the 
macroscopic and microscopic part of electronic density 
response and the perturbing potential, via Eqs. (45), 
(39), (41), and (47); then solve Eq. (46) for A E V> and 
repeat. However, for computational purposes, it is more 
convenient and simple to keep the value of the screened 
electric field, E, fixed and let only the microscopic compo- 
nents of the potential vary during the the iterative proce- 
dure, that is reached simply iterating Eq. (46), (41) and 
(47). The macroscopic polarization is then calculated 
only at the end when self-consistency is achieved using 
Eq. (45). Physically, this is equivalent to calculating the 
polarization response to a given screened electric field, E, 
instead of the response to the bare electric field E . 

The electronic contribution to the dielectric tensor, 
e£f , for the general (low-symmetry) case, can be derived 
from simple electrostatics. Using Eq. (39) and the defi- 
nition of e£f , one has: 

E 0a = (E a +47rP a )=^eSfE /3 . (48) 

P 

Using Eq. (45) to calculate the polarization induced in 
the a direction when a field is applied in the (3 direction, 
one finally obtains: 

N/2 



F a/3_* 167Fe V^aiAE, 



H n=l 



(49) 



3. Relation to the variational principle 

The KS equations are the Euler equations which solve 
the Hohcnbcrg-Kohn variational principle. The equa- 
tions of DFPT introduced in Sec. II. C can be seen as 
a set of equations which solve approximately the varia- 
tional principle when the external potential is perturbed. 
Alternatively, these equations can be considered as those 
which minimize exactly an approximate energy functional 
(Gonzc et al. 1992, Gonze 1995a, Gonze 1997). To see 
this point in some detail, let us consider the energy func- 
tional as depending explicitly on the set of KS orbitals, 
ip = {ip n } (assumed to be real) and parametrically on 
the external potential, V(r): 

ra=l 

-J V (r)n(r)dr + ^ j 



+ 



drdr' 



r - r'| 

+E xc [n(T)]. (50) 



The functional derivative of the above functional with 
respect to ip n is: 



5E 



2H SCF ip n (v). 



Stp n (r) 

The Euler equations thus read: 



N/2 

Hscf\^Pti) — ^ ' A m „ 
m—1 



\lpm), 



(51) 



(52) 



where the A's are a set of Lagrange multipliers intro- 
duced so as to enforce the ortho-normality of the ip's. 
Eqs. (52) are invariant with respect to a unitary transfor- 
mation within the manifold of the ip's, so that the usual 
KS equation, Eq. (15), is recovered in the representation 
which diagonalizes the A matrix. 

Let us now indicate by ip( > the solutions of the KS 
equations corresponding to a particular choice of the 
external potential, V^(r) (the unperturbed potential), 
and let us indicate by AV and Aip the differences be- 
tween the actual potential and orbitals and their un- 
perturbed values. The energy functional, Eq. (50), can 
be equally seen as depending on Aip and AV: E = 
E[{tpW + Aip};V^ + AF]. We now consider the ap- 
proximate functional, E^ , which is obtained from E by 
truncating its Taylor expansion in terms of Aip and AU 
to second order: 



£ (2) [{AiP}; AV] = E[{iP^}- V<°>] - 



SE 
SV(r) 



AV(r)dr 



8E 



N/2 

§7 SMr) 



Aip n (r)dr 
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N/2 



■It 



s 2 e 



2^ J 5Mr)SV(r') 

N/2 



-A^ n (r)AV(r')drdr' 



+ 



\Tj s^f^ WMWm 



where we have omitted the second derivative with respect 
to V because E is linear in V. Besides Eq. (51), the 
required functional derivatives are: 



SE 



8V(r) 



S 2 E 



SMr)SV(r') 
5 2 E 



5tp n (r)5tp m (r') 



i/.=i/ ,<0> 



(54) 
(55) 



n(°)(r) 

2^(r)<5(r-r') 
+4X(r,r')^ 0) (r)^(r'), (56) 



where: 



K(r,r') 



S 2 E XC 



r — i" 



5n(r)Sn(r') ' 



(57) 



and ft4°J F (r, r') is the kernel of the unperturbed KS 
Hamiltonian. The energy functional, Eq. (53), must be 
minimized under the constraints that the resulting solu- 
tions lead to an orthonormal set of occupied states: 



This leads to the Euler equations: 



(58) 



N/2 



m — 1 



N/2 



-AV SCF \^) + J2 ( A nm - e n S nm )\^), (59) 



m— 1 



where the A's are a set of Lagrange multipliers introduced 
so as to enforce the constraints, Eq. (58), and 

AV SCF (r) = AV(r) 

N/2 

+2J2 dr'K(r,r')^(r')AMr') (60) 

71=1"' 

is the first-order variation of the self-consistent potential. 
We now project both sides of Eq. (59) onto ip^ . Taking 
into account that by Eq. (58) (A^ n |^ 0) ) = C(A 2 ), to 
linear order in A one obtains: 



av scf \^). 



(61) 



To linear order in A, Eq. (59) can thus be cast into the 
form: 



(Hi% - e n )|AV„> = -AV SCF \^) 

N/2 

+ £iva<vaA^ivi o) > 



rn — 1 



-PcAVscfI^), 



(62) 



which is essentially the same as Eq. (30). The set of 
perturbed orbitals, {Aip}, is the solution of N coupled 
linear systems whose dimension is the size of the basis 
set, M (Eq. (62) or Eq. (30), where the coupling comes 
from the dependence of the self-consistent potential on 
all the orbitals). Alternatively, it can be seen as a huge 
linear system of dimension MxN/2, which is obtained by 
inserting the expression for AV^ SC f, Eq. (60), into Eqs. 
(62) or (30). The latter is naturally derived from the 
minimization of the functional E^ 2 \ Eq. (53), which is 
quadratic. 

This variational approach shows that the error on the 
functional to be minimized, Eq. (53), is proportional to 
the square of the error on the minimization variables, 
Aip. This fact can be exploited in the calculation of the 
second-order mixed derivatives, Eq. (22). It can be shown 
that a variational expression can be constructed for 
mixed derivatives as well (Gonze et al. 1992, Gonze 1997). 



4. Metals 

Density-functional perturbation theory, as presented 
above, is directly applicable to metals, provided the (elec- 
tronic) temperature vanishes, so that a clear cut sepa- 
ration between occupied and empty states is possible. 
In this case, however, the number of k points needed 
to correctly represent the effect of the Fermi surface 
would be very large. Practical implementations of DFPT 
to metallic systems have been discussed by Quong and 
Klein (1992) and by de Gironcoli (1995), in the pseudo- 
potential formalism, and by Savrasov (1992), in the Lin- 
earized Muffin-Tin Orbital framework. In the follow- 
ing, we will closely follow the formulation of de Giron- 
coli (1995) which is based on the smearing technique for 
dealing with Fermi-surface effects. 

In the smearing approach each KS energy level is 
broadened by a smearing function: 



(63) 



where S(x) is any function that integrates to 1 — so that 
<5 CT (e) tends to the Dirac S function in the limit of van- 
ishing smearing width, a. Many kinds of smearing 
functions, 6(x), can be used: Fermi-Dirac broadening, 3 



In this case the broadening function is the derivative of the 
Fermi-Dirac distribution function: 8(x) = | (1 + cosh(x)) -1 . 
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Lorentzian, Gaussian (Fu and Ho 1983), Gaussian com- 
bined with polynomials (Methfessel and Paxton 1989), 
or cold smearing (Marzari et al. 1999) functions, to recall 
just a few of those used in the literature. While the choice 
of a given smearing function is to some extent a matter of 
personal taste and computational convenience, the spe- 
cific choice of the Fermi-Dirac broadening allows one to 
explicitly account for the effects of a finite temperature 
(T = er/fcs), when needed. The (local) density of states 
resulting from the broadened energy levels will be the 
original density of states, convoluted with the smearing 
function: 

n(r,e)=]TW^W(r)| 2 , (64) 

where the sum refers both to the discrete k-vector index 
and to band and spin indices, for all bands. From this 
basic quantity the electron density follows: 

n(r)= J e \(r,e)de = ^9(^—^) |^„(r)| 2 , (65) 

where 9(x) = S(y)dy is a smooth approximation to 
the step function, and the Fermi energy is determined by 
the normalization to the total number of electrons, 

N = f F n(e)de = ^ (tZJllA . ( 66 ) 
J-oo n \ a ) 

The advantage of this procedure is that after convolution 
the (modified) local density of states can be computed ac- 
curately on a discrete grid of points in the BZ, provided 
that the average separation between neighboring eigen- 
values is small with respect to the broadening width, a. 

Given the definition of the local density of states, 
Eq. (64), the consistent way to define the auxiliary KS 
kinetic-energy functional — or its analogue in the finite 
temperature theory (Mermin 1965) — is through the Leg- 
endre transform of the single-particle energy integral: 




where Q\ (x) — yS(y)dy. Note that when a finite elec- 
tronic temperature is considered the KS auxiliary func- 
tional, T s [n], contains both the kinetic energy and the 
entropy contribution to the electronic free energy of the 
independent particle system. These two contributions 
appear separately in the last expression for T s [n] in Eq. 
(67) where it can be verified that for Fermi-Dirac broad- 
cning h{x) = /log(/) + (1 - /)log(l - /), with / = 6{x) 



as required. With the definitions above, for any kind of 
smearing function, the usual KS equations follow from 
the minimization of the total energy. The price to be 
paid for the computational simplicity of the smearing ap- 
proach is that the computed total energy depends on the 
chosen broadening width and results for finite broaden- 
ing widths must be corrected for, unless the shape of 
the smearing function is such as to reduce this depen- 
dence within acceptable values. Various discussions of 
this issue can be found in the literature (Methfessel and 
Paxton 1989, Marzari et al. 1999, de Gironcoli 1995, De 
Vita 1992). 

Forces and other first-order energy derivatives can be 
computed in the usual way from Eqs. (8) and (21) which 
only require the knowledge of the unperturbed electronic 
density. Similarly second-order derivatives are computed 
from Eq. (22) where also the first order variation of the 
density is needed. Direct variation of Eq. (65) provides 
the required expression for An(r): 

An(r) = ^> [C«AV„(r) + c.c] 

n 

+ ^|^„(r)| 2 ^„(Ae F - Ae„), (68) 

n 

where we have defined 9 n , m = 9((e n — e m )/a) and 

S n ,m = (l/ <7 )^(( £ n — fm)/^- The last term in Eq. 
(68) accounts for possible changes in occupation num- 
bers induced by shifts in the single particle energies 
(Ae„ = (ip n \AV SCF \ip n )), as well as in the Fermi energy 
of the system. Whether this term is present or not de- 
pends on the thermodynamic ensemble used: it is absent 
if the chemical potential is kept fixed, whereas it might 
be present if it is the number of electrons which is fixed. 
Even in this last case the Fermi energy is unaffected by 
the perturbation to linear order, unless the perturbation 
is lattice-periodic (monochromatic with q = 0). Let us 
neglect this term for the time being and come back to it 
at the end of the section. 

Substituting in Eq. (68) the definition for A-0„(r) from 
standard perturbation theory, Eq. (28), and exploiting 
the symmetry between the two contributions in square 
brackets, one obtains the following well known expression 
for the first-order variation of the electronic density in a 
metal: 

a / \ @F,n ~ 0F,m 

An(r) = ^ — — 

x V;(r)Vm(r)(Vm|Ay scF |^„), (69) 

where the term in Ae„ in Eq. (68) has become the 
n = m term in the above sum and the incremental ra- 
tio (0F,n — 0F,m)/(en — £m) must be substituted with its 
limit, — Sf,u, whenever e m — > e n . This limit is always fi- 
nite for any finite broadening line- width, or temperature, 
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and this expression is therefore numerically stable even 
in the presence of vanishingly small virtual excitation en- 
ergies. 

In order to avoid the double sum over occupied and 
unoccupied states we use the relation 6{x) + 0(—x) = 1 
and the symmetry between i and j to get 



An(r) = 2]T 



x V;(r)Vm(r)(Vm|AV SCF |^ n >, (70) 



where the sum over the first index can be limited to the 
states that have non negligible occupation. This expres- 
sion can be further simplified, avoiding the explicit sum 
over the second index, by re- writing it as: 



An(r) = 2^C«A^(r), 



(71) 



where the A^„'s satisfy the equation 



Q 



Aip n 



'F,n 



Pn 



Ay sCF |V„), (72) 



with 



(73) 



In the above equations the a^'s are chosen in such a way 
that the Q operator makes the linear system, Eq. (72), 
nonsingular for all non vanishing Aip n . A possible simple 
choice is 



a k = max(e F + A - e fe , 0), 



(74) 



with A w 3 -j- 4(7. Another, even simpler, choice is to 
set ak equal to the occupied bandwidth plus, say, 3<r for 
all partially occupied states and zero when the state is 
totally unoccupied. It can be easily verified that since 
vanishes when ipk is unoccupied then also (3 n;m vanishes 
when any of its indices refers to an unoccupied state. 
Therefore the Q and P operators involve only the small 
number of partially filled bands and the first-order varia- 
tion of the wave functions and of the charge density can 
be computed avoiding any explicit reference to the un- 
occupied states, much in the same way as for insulating 
materials. In fact, if the above scheme is applied to an 
insulator using a smearing width much smaller than its 
fundamental band gap, all the metallic equations, Eqs. 
(71-73), reduce numerically to their insulating analogues, 
Eqs. (23) and (30). 

The expression for the charge density linearly induced 
by a given perturbation, Eqs. (69) and (70), involves an 



energy denominator which vanishes for metals. In one 
dimension, this vanishing denominator gives rise to a di- 
vergence in the screening of perturbations whose wave- 
vector is twice the Fermi momentum, 2kp. This di- 
vergence is smeared in two dimensions and suppressed 
in three dimensions by volume effects. However, if the 
topology of the Fermi surface is such that two finite por- 
tions of it are parallel and connected by a wave-vector 
which for convenience we will name 2fcp {nesting), the 
screening to perturbations of wave-vector 2kp will di- 
verge even in three dimension. This is the physical mech- 
anism giving rise to Kohn anomalies in the vibrational 
spectra of certain metals. The sampling of the BZ neces- 
sary to evaluate the sum over (partially) occupied states 
in Eqs. (69) and (70) is in ordinary cases similar to that 
needed to calculate the unperturbed charge-density dis- 
tribution. Near a Kohn anomaly, however, a fine sam- 
pling of the Fermi surface is necessary and the number 
of needed points in the BZ is correspondingly larger. 

Periodic (q = 0) perturbations may induce a shift of 
the Fermi energy. In this case Eq. (71) must be modified 
as: 

An(r) = 2]TC(r)A^(r)+n(r,e F )Ae F , (75) 

n 

all other technical details remaining unchanged. In 
order to find out the appropriate value of the Fermi 
energy shift let us examine the perturbation in the 
q — > limit. Let us consider the Fourier transform 
of the self-consistent perturbing potential Ay sCF (q) = 
1/V J AV r SCF (r)e _iq ' r dr. Its macroscopic (q w 0) com- 
ponent reads: 

AV SCF {q) = AV(q) + — An(q) + ^An(q), (76) 
q z an 

where the last term is the exchange-correlation contribu- 
tion, AV(q) = — ^§-An ext (q) is the macroscopic elec- 
trostatic component of the external perturbing potential, 
and An ext (q) is finite in the q — > limit. On the other 
hand the macroscopic component of the density response 
is 

An(q) = -n(e F )AV SCF (q) + An lf (q), (77) 

where An 1 * is the density response to the non- 
macroscopic (local fields) components of the self- 
consistent potential. As An 1 * and the DOS at the Fermi 
energy, n(e F ), are both finite, AV SCF and An must not 
diverge when q — > otherwise Eqs. (76) and (77) could 
not be satisfied at the same time. This implies, from 
Eq. (76), macroscopic charge neutrality for the perturbed 
system, that is An ext (q) = An(q) + 0(q 2 ) for q w 0; a 
condition that in turn implies in Eq. (77) 



AF SCF (q) 



An e!Ef (q)- An^(q) 
n(e F ) 



+ 0(q 2 ). (78) 
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The lattice periodic (q = 0) result can be obtained by 
taking the q — > limit of the above equations. In 
this case however, it is costumary to set arbitrarily to 
zero the macroscopic electrostatic component of the self- 
consistent potential and the charge neutrality condition 
is thus enforced by a compensating shift in the Fermi 
energy, equal and opposite to the above result: 

An eit (q = 0)+/n(r,e f )AV(r)* 

Ae_p = — . (79) 

n(e F ) 



D. Phonons 

1. Vibrational states in crystalline solids 

In crystalline solids, the nuclear positions appearing in 
the definition of the IFC's, Eq. (5), are labelled by an 
index, /, which indicates the unit cell to which a given 
atom belongs, I, and the position of the atom within that 
unit cell, s: I = {l,s}. The position of the 7-th atom is 
thus: 



Rj = Ri+t 8 + 11,(0, 



(80) 



where R; is the position of the Z-th unit cell in the Bravais 
lattice, t s is the equilibrium position of the s-th atom in 
the unit cell, and u s (l) indicates the deviation from equi- 
librium of the nuclear position. Because of translational 
invariance, the matrix of the IFC's, Eq. (10), depends on 
I and m only through the difference R = R; R m : 



CS"(i,m) = 



d 2 E 



C 5 f(R,-R ro ), (81) 



duf(l)du^(m) 



where the greek super-scripts indicate Cartesian compo- 
nents. The Fourier transform of Cff (R) with respect to 
R, C"f(q), can be seen as the second derivative of the 
BO energy surface with respect to the amplitude of a 
lattice distortion of definite wave-vector: 



C s f(q) 



V" g-iq-R-rr"/ 3 

R 



c s f(R) 



d 2 E 



N c du^(q)du f ^ q y 



(82) 



where N c is the number of unit cells in the crystal, and 
the vector u s (q) is defined by the distortion pattern: 



R/[u s (q)] = R/ +t s + u s (q)e lqRi . 



(83) 



Phonon frequencies, w(q), are solutions of the secular 
equation: 



dct 



1 



C s f(q)-^(q) 



0. 



(84) 



Translational invariance can be alternatively stated in 
this context by saying that a lattice distortion of wave- 
vector q does not induce a force response in the crystal at 
wave-vector q' ^ q, in agreement with the analysis car- 
ried out in Sec. II.C.l. Because of this property, IFC's 
are most easily calculated in reciprocal space and, when 
they are needed in direct space, they can be readily ob- 
tained by Fourier transform (see Sec. II. D. 3). 

The reciprocal-space expression for the matrix of 
IFC's, Eq. (10), is the sum of an electronic and of an 
ionic contribution: 



(q) = ei C,T(q) + 4 °"C s T(q), 



_el s-ia.fi 



(85) 



where: 



st (q) = w 



/G 



9n(r) \* dV ion (r) 



dr 



n(r) 



d 2 V lon {r) 
0u;°(q)0u£(q) 



-dr 



and 



V ion (r) 



v s (r - R; - T s - u s (l)), 



(86) 



(87) 



v s being the ionic (pscudo-) potential corresponding to 
the s-th atomic species. All the derivatives must be cal- 
culated for u s (q) = 0. The ionic contribution comes from 
the ion-ion interaction energy (the last term of Eq. (10)) 
and does not depend onJ;he electronic structure. The 
explicit expression of 4 °"C"/(q) for periodic systems is 
given in the appendix. 

Using Eqs. (83) and (87), the derivatives of the poten- 
tial appearing in Eq. (86) read: 

dV ion {r) = _ flw.fr-Rt -T,) c<q . Rl 



0u?(q) 



dr 



while the corresponding derivative of the electron charge- 
density distribution is given by Eqs. (33) and (35). 



2. Long wave-length vibrations in polar materials 

In polar semiconductors and insulators, the long-range 
character of the Coulomb forces gives rise to macroscopic 
electric fields for longitudinal optic (LO) phonons in the 
long wave-length limit. At any finite wave-length, polar 
semiconductors are dealt with in the same way as non po- 
lar ones. In the long wave-length limit, however, phonons 
are coupled to macroscopic electric fields which must be 
treated with some care because the corresponding elec- 
tronic potential, V^(r) = eE • r, is not lattice-periodic 
(see Sec. II. C. 2). A physically transparent picture of the 
coupling between zone-center phonons and macroscopic 
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electric fields is provided by the Huang's phcnomcnolog- 
ical model (Born and Huang 1954), which we briefly dis- 
cuss in the case of a cubic (or tetrahedral) lattice with 
two atoms per unit cell. The most general quadratic ex- 
pression of the energy as a function of the phonon optic 
coordinates, u, and the electrical degrees of freedom (i.e. 
the field itself, E), is: 

E(u,E) = \mujIu 2 - ^- eoo E 2 

where M is the nuclear reduced mass, f2 is the volume of 
the unit cell, eoo the electronic dielectric constant of the 
crystal (i.e. the static dielectric constant with clamped 
nuclei, u = 0), and the coupling between u and E, Z*, 
is known as the Born effective charge of the ions (see 
e.g., Bottger (1983), Sect. 1.5). The variables which are 
conjugate to u and E are the force F acting on the ions 
and the electrical induction, D: 



Z*u ■ E, 



(89) 



-Mu 2 u + Z 
Aw 



dE 
^ ~du = 
AwdE__ 

In the absence of free external charges, Maxwell equa- 
tions give: 



Z*u + e OQ E. 



(90) 
(91) 



rot E <~ iq x E = 
div D ~ iq • D = 0. 



(92) 
(93) 



For transverse modes 
and Eq. (90) F T = - 
is therefore uot = ^o- 



(E 1 q), Eq. (92) gives E T = 0, 
MwqU: the transverse frequency 
For longitudinal modes (E || q), 



Eq. 
Eq. 



(93) gives D L 
(90) gives F L 



and Eq. 

= - (MloI - 



(91) gives 



Sle = 



u: 



4ttZ* 

•-l - --n^"> 

the longitudi- 



4ttZ* 



These 



nal frequency is therefore u>l = \J ^5 ~<~ He — m- 
results, which are exact in the case of cubic and tetra- 
hedral systems, can be easily generalized to crystals of 
arbitrary symmetry (Bottger 1983). 

The first-principles calculation of e M and Z* proceeds 
from Eqs. (90) and (91). Let us start for instance from 
Eq. (91) which — expressed in terms of the macroscopic 
electric polarization of the medium and generalized to 
the case of many atoms per cell — reads: 



i 



Aw 



(94) 



In the general, low-symmetry case, Eq. (94) must be 
read as a tensor equation stating that the Born effective- 
charge tensor of the s-th ion is the partial derivative of 
the macroscopic polarization with respect to a periodic 
displacement of all the ions of the s-species at zero macro- 
scopic electric field: 



z* a / = n 



du l3 s {q = 0) 



(95) 



while the electronic dielectric-constant tensor is the 
derivative of the polarization with respect to the macro- 
scopic electric field at clamped nuclei: 



e a J = 5 a p + Aw 



dP 



dE, 



(96) 



u 3 (q=0)=0 



In the long wave-length limit, the matrix of the force 
constants can be split into the sum of an analytic and 
a non analytic contribution (Born and Huang 1954, 
Cochran and Cowley 1962): 



s-ia.fi _an /~ia(3 . na stotf) 



(97) 



where the analytic part, a "C, is the matrix obtained from 
the response to a zone-center phonon, calculated at zero 
macroscopic electric field. The non analytic part has the 
general form (Cochran and Cowley 1962): 



Aw (q-Z%) a (q-Z* t ) /3 



q • C°° • q 



(98) 



Eq. (98) shows that all the information necessary to deal 
with the non analytic part of the dynamical matrix is 
contained in the macroscopic dielectric constant of the 
system, and in the Born effective charges Z* , whereas 
the analytic contribution can be calculated just ignoring 
any macroscopic polarization associated to the phonon. 
All these quantities can be easily obtained within DFPT 
(Giannozzi et al. 1991). 

It is worth mentioning that effective charges can be cal- 
culated using an approach to the electrostatics of quan- 
tum dielectrics based on topological concepts, the Berry 's 
phase approach to macroscopic polarization (King-Smith 
and Vanderbilt 1993, Resta 1994). When used at the 
same level of accuracy, the linear-response and Berry's 
phase approaches yield the same results within numeri- 
cal uncertainties. 



3. Inter-atomic force constants 

The considerations developed so far in the present sec- 
tion allow in principle to calculate the vibrational fre- 
quencies at any (finite or infinite) phonon wave-length. 
Phonon frequencies are usually rather smooth functions 
of the wave-vector, so that suitable interpolation tech- 
niques can be used when complete dispersions arc needed. 
Simple concepts from (discrete) Fourier analysis show 
that the smoother the phonon dispersions (i.e. the 
smoother the matrix elements of C as functions of q), 
the shorter is the range of real-space IFC's: 



E=o 



C s f(R) = ^^ e ^ R C s f(q), 



(99) 
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i.e. the smaller the number of their non- vanishing val- 
ues (to any given accuracy). Real-space IFC's can thus 
be readily obtained by Fourier analyzing a set of force- 
constant matrices calculated and tabulated over a uni- 
form grid of points in reciprocal space. The most efficient 
way of calculating all these Fourier transforms numeri- 
cally is the fast Fourier transform (FFT) technique (see 
e.g. Press ct al. (1989)). Once real-space IFC's have 
been thus obtained, dynamical matrices in reciprocal- 
space (and, hence, vibrational frequencies) can be ob- 
tained at any wave- vector (not necessarily contained in 
the original grid) by FFT. The shorter the range of real- 
space force constants, the coarser will be the reciprocal- 
space grid needed for such Fourier interpolation. In prac- 
tice, the size of the reciprocal-space grid will be assessed a 
posteriori by verifying that it yields vanishing real-space 
constants (to within a given accuracy) beyond some cut- 
off radius. A simple rule of thumb is to include in the 
FFT grid enough points in the BZ so as to reach neigh- 
bor interactions extending up to 2-3 bond- lengths, and 
check the accuracy of the interpolation against the full 
calculation on some points not included in the grid. 

The above considerations apply as such to metals, 
away from Kohn anomalies, and to non-polar insulators. 
The presence of Kohn anomalies in metals is associated 
to long-range IFC's propagating along the direction of 
the wave- vector of the anomaly. Catching the details of 
the anomaly with a regular grid of wave- vectors in the BZ 
would be very impractical. In these cases, once the posi- 
tion of the anomaly has been located, it is much simpler 
to refine the grid locally just around the anomaly. 

In polar materials, instead, real-space IFC's are long- 
ranged in all directions, as a consequence of the long- 
range dipole-dipolc interaction between ionic effective 
charges. For this reason, Fourier interpolation would be 
inefficient in this case. The dipole-dipole interaction is 
precisely the physical origin of the non-analytic behavior 
of the reciprocal-space dynamical matrices in the long 
wave-length limit, whose form is however known exactly 
in terms of the ionic effective charges (Eq. (98)). Eq. (98) 
expresses the long wave-length limit of the reciprocal- 
space force constants of any system whose atoms carry a 
charge equal to Z* s . If the force constants of a system 
of point charges Z* s is subtracted from those calculated 
for the physical system under consideration, the resulting 
difference will be analytic in the long wave-length limit, 
and its Fourier transform short ranged. For polar ma- 
terials, Fourier interpolation is thus efficiently applicable 
to the analytic contribution to the reciprocal-space force 
constants, whereas the full non-analytic behavior can be 
easily restored by adding the force constants of a suitable 
point-charge model (Giannozzi ct al. 1991). A descrip- 
tion of the technical details necessary to implement the 
Fourier interpolation of reciprocal-space force constants 
in the general case of materials with anisotropic effective 



charges can be found in Gonze and Lee (1997). 

E. Homogeneous deformations 

1. Elastic properties 

Elastic constants can be viewed as force constants as- 
sociated to homogeneous strains, i.e. to macroscopic dis- 
tortions of the crystal. In any finite system, there is no 
conceptual difference between a strain and a microscopic 
distortion, and linear-response techniques are straightfor- 
wardly applicable in both cases. In an infinite system, on 
the contrary, one cannot apply directly linear-response 
techniques, because a homogeneous strain changes the 
boundary conditions of the Hamiltonian. The use of per- 
turbation theory requires instead the existence of a com- 
mon basis set for the perturbed and unperturbed sys- 
tems. In order to use perturbation theory for homoge- 
neous deformations, it has been suggested that one can 
introduce an intermediate fictitious Hamiltonian which 
is related to the unperturbed one by an unitary transfor- 
mation, and obeys the same boundary conditions as the 
strained Hamiltonian (Baroni ct al. 1987b). 

Let us consider for simplicity a system under an 
isotropic strain (dilatation) of amplitude a, {R} — > 
{aR}, whose corresponding clastic constant is the bulk 
modulus. The KS Hamiltonian, H a , for the strained 
crystal is given by 

h 2 d 2 

SCF 2mdr 2 lon[) 

+ e 2 J yf^7]dv' + v xc (n a (r)), (100) 

where the electron- ion potential, V^ n , is 

V? on {v)=Y J V s {v- a K l - a T s ). (101) 

Is 

The intermediate (fictitious) strained Hamiltonian, 
H a , is obtained from the unperturbed one, H, through 
a scale transformation: 

H a (r, d/dv) = H SCF {v/a, ad/Or). (102) 

H a obeys the same boundary conditions as the physical 
strained Hamiltonian H a , and hence perturbation theory 
can be used to calculate the relative energy difference. At 
the same time H a and H differ by a unitary transforma- 
tion and their spectra are trivially related: 

C(r) = a-^ n (v/a) (103) 
n a (r) = a^n{v/a). 
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The energy change due to a strain can thus be com- 
puted in two steps: first one calculates the energy dif- 
ference between the unperturbed crystal and the ficti- 
tious strained one described by H a ; one then computes 
the energy difference between the latter and the physi- 
cal strained system using perturbation theory. The first 
step is trivial. The second step is less so, because H a is 
not a proper Kohn-Sham Hamiltonian: the Hartree and 
XC terms in H a are not the Hartree and XC potentials 
generated by n a . One can write H a as a genuine Kohn- 
Sham Hamiltonian, provided that the definition of the 
external potential is changed: 



where 



h 2 a 2 d 2 



2m dr 2 +V ^^ 



h a (r') 



j-dr' + v xc {h a (r)), (104) 



K a xt (r)=5> s (r/a-R ; -T s 



1-1 

a 



rc(r') 



r/a - 



:dr' 



+ v xc (a 3 n(r/a)) + v xc (n(r / a)) . (105) 

The energy difference between the fictitious and the real 
strained systems, E a — E a , can now be calculated by 
perturbation. The details are given in Baroni et al. 
(1987b). This method can be straightforwardly extended 
to generic elastic constants. The algebra involved is how- 
ever quite heavy. 

It must be remarked that what is calculated in this 
way is only a bare elastic constants: in Eq. (101) the co- 
ordinates of different atoms within a same unit cell are 
assumed to undergo the same homogeneous scale trans- 
formation as the positions of different unit cells. This 
is not true in general, and atoms rearrange themselves 
within each unit cell, so as to minimize the total energy 
as a function of the applied strain. To see how the elas- 
tic constants are affected by the internal relaxation of the 
atomic positions, let us write the more general second- 
order expression of the crystal energy per unit cell as a 
second-order polynomial of the macroscopic (strain) and 
microscopic (phonon) modes: 

n 



l 

2 



macroscopic strain and atomic displacements (the ma- 
trix of the so-called internal- strain parameters) . Crystal 
symmetry determines the number of independent non- 
vanishing terms in A , C, and £. If we allow the atoms 
to relax for a given strain state, minimization of the en- 
ergy with respect to the w's yields 



0,-yS Ca0 ^7<5i 



(107) 



a/3, 7(5 



where 



Aa/3,7,5 — ~ ^ Csfj.a/3 (C ^gtCtv-yS- (108) 

Both the force constants and the coupling between 
macroscopic strain and atomic displacements can be 
readily calculated within DFPT. 



2. Piezo-electric properties 

The piezo-electric constants form a third-order tensor, 
7a, 75, defined as the derivative of the macroscopic elec- 
tric polarization with respect to a homogeneous strain, 
at vanishing macroscopic field. This quantity — which has 
been demonstrated to be independent of surface effects 
(Martin 1972), i.e. independent on surface termination — 
could in principle be evaluated as the electric-polarization 
response to a given applied strain. Alternatively, and 
somewhat more conveniently, it can also be calculated 
as the stress linearly induced by an electric field at zero 
strain. The latter definition was used by de Gironcoli 
et al. (1989) to calculate the piezo-electric tensor in III- 
V semiconductor compounds. The two definitions are 
equivalent and can be deduced expanding the energy of 
the system to second order in the macroscopic perturba- 
tions (strain, e a p, and electric field, E a ): 

n 



Si X ^ 

E = — A a/ 3 i7 a e a p e 7< 5 



a/3, 7(5 



n 



-Sl^ 7ai75 E a e jS - — ^ef E a E p . (109) 



a, 7(5 



a,/3 



The link with the microscopic description is provided, 
by expanding the crystal energy to second-order in the 
external (strain and electric field) and internal (phonon) 
degrees of freedom: 

ft 



+ C "f <u^ + Qj2 C~7« < e 7* (106) E = - ]T \% M e a p e^ s ~nY, 7".7« E « e 7* 



sa,t/3 



sajS 



where Si is the unit cell volume, A is the (bare) elastic 
constant matrix, e the strain tensor, C the zone-center 
reciprocal-space matrix of force constants, u atomic dis- 
placements in the unit cell, Q the coupling between 



a/3, 7(5 

~~ e °° tct h P + 9 ° s * u * Ut 



a,7<5 



+ SU7J«?e7!-flj;^' ( 110 ) 

sajS safd 
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where A , 7 , and are the purely electronic (i.e. 
clamped-ion) contributions to the elastic, piezoelectric, 
and dielectric, constants respectively, and the coupling 
between atomic displacements and macroscopic variables 
(electric field and strain), are expressed by the effective 
charges, Z* , and internal strain parameters, (. Once the 
macroscopic variables, E and e, are fixed, The equilib- 
rium values of the internal degrees of freedom are given 
by the condition that the derivatives of Eq. (110) with 
respect to the atomic displacements vanish. When these 
equilibrium atomic positions as functions of the macro- 
scopic electric field and strain are inserted in (110), the 
resulting expression defines the total piezoelectric con- 
stants as: 



E 



^ s 1° 1st SytwyS- 



(111) 



The two resulting contributions to the piezoelectric con- 
stants are often of opposite sign and close in absolute 
value, so that a well converged calculation is needed in 
order to extract a reliable value for their sum (dc Giron- 
coli et al. 1989). 

The problem of a proper definition of piezoelectric 
properties in a crystal displaying a spontaneous macro- 
scopic polarization has been raised recently by Saghi- 
Szabo et al. (1998) and further discussed by Vanderbilt 
(2000). 



F. Higher-order responses 

1. The 2n+l theorem 

In Sec. II. C we have seen that the knowledge of the 
first-order derivatives of the wave-functions is enough to 
calculate the second-order derivatives of the total energy. 
This a special case of a very general theorem, known as 
the 2n+ 1 theorem, which states that the knowledge of the 
derivatives of the wave-functions up to order n allows the 
calculation of the derivatives of the energy up to order 
2n + 1. This theorem, well known in quantum mechan- 
ics since many years, is a consequence of the variational 
principle and it is valid also in DFT. In this context, 
its usefulness derives from the fact that the third-order 
derivatives of the total energy can be obtained from the 
first-order derivatives of the wave- functions. This opens 
the possibility to study phenomena which depends upon 
third order anharmonic terms in the energy expansion — 
such as phonon line widths, Raman scattering cross sec- 
tions, or nonlinear optical responses — with a computa- 
tional effort of the same order as for harmonic properties 
(because the time-consuming step is the calculation of 
the first-order derivatives of the wave-functions). 

Several proofs of the 2n + 1 theorem can be found in 
the literature. In a DFT framework this theorem was 



first proved by Gonze and Vigneron (1989). Explicit 
expressions of the energy derivatives up to fourth order 
have been worked out by Gonze (1995b). The scope of 
this theorem is much more general, as it concerns the 
properties of the extrema of any functional depending 
upon some parameters (Epstein 1974). Let E[tl>, A] be 
a generic functional of ip which, for A = 0, has an ex- 
tremum at V (0) : 5E[il> ia) ,0}/5i> = 0. The position of the 
extremum will depend on the value of the parameter A: 
V>(A) = tp^ + Aip(\). The value of the functional at the 
extremum will be: 



E min (\)=E[4> {a) +Atl>(\),\], 



(112) 



where A^>(A) is determined, for any given value of A, by 
the extremum condition: 



dE[xP {0) + AV>, A] 
dAtj) 



= 0. 



(113) 



When A is small, both Aip(X) and E min (X) will be well 
approximated by their Taylor expansion in powers of A: 



1=1 



(114) 



1=1 



iWA) = £^%^ = E s( ° A - ( 115 ) 

1=0 ' 1=0 

The 2n + 1 theorem states that the knowledge of At/}^ 
up to order n is enough to determine E^ up to order 
2n + 1. To demonstrate this, it is convenient to first 
expand E[tp^ + Alp, A] into a Taylor series treating Ai{> 
and A as independent variables: 

f,(0) 



£[V> W +AV>,A] 



where we use the notation: 



^|(A^) fc 



(116) 



(117) 



Variation of Eq. (116) with respect to Alp leads to the 
extremum condition: 

dE[^ 0) +AV>,A] 
dAtp 



= 0. 



(118) 



Introducing Eq. (114) into Eq. (118) allows to formally 
expand f as a power series of A only. The resulting ex- 
pression for f must vanish identically order by order in 
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A. By equating the coefficients of A™ in the Taylor ex- 
pansion to zero, Eq. (118) yields an infinite number of 
equations, f(") = 0, which determine Atp^ n \ 

Introducing Eq. (114) into Eq. (116) allows to formally 
expand E as a power series of A only and therefore to cal- 
culate the explicit form of E®. Since a term quadratic 
in Aip^X 1 is of order A 2 ', E^ 2n+1 ^ could depend only lin- 
early upon A^/^A', if I > n. Showing that for I > n the 
coefficient of A^ l) X l in E (2n+1 1 is f( 2 "+!-0 we prove the 
2n + 1 theorem, since this force is zero by the minimum 
condition. We start by extracting Ai()^X l from the prod- 
uct (Aij)) appearing in Eq. (116), using the relationship: 



(A^) te = (AV> - Aip {l) X l + Aip (l> X l ) k = (119) 
kA^X^A^f- 1 + {Alp - Ai> {l) X l ) k + 0(X 2n+1 ), 

valid for I > n. The only term which is linear in A^i^A' 
is the first term of the r.h.s. of Eq. (119). Inserting this 
term in Eq. (116) and recalling the definition of f, Eq. 
(118), we can write E^+V as: 

E (2n+1) = A ^(2«+l) A 2n+l f (0) + ^ 

+ A^ l) X l f( 2n+1 - l > + . . . 

+ A^ (n+1) A (n+1) f (n) + p( 2 " +1 )(V> (1) , . . .,V (n) ), 

where p( 2 ™+ 1 ) is a polynomial of degree 2n + 1. From 
the condition fW = for every i, we get: 



E (2n+1) = p(2n+l)(^(l) } ... i ^(»)) i 



(121) 



that proves the 2n + 1 theorem. 

In order to apply this theorem to DFT we can take as 
%j) a vector whose elements are the coefficients of all the 
occupied wave- functions {ipi} on a given basis and A as a 
parameter measuring the magnitude of the perturbation. 
The orthogonality constraint can be dealt with as shown 
for instance by Mauri et al. (1993), writing the total en- 
ergy functional of DFT for non ortho- normalized orbitals. 
The demonstration presented here, applied to this func- 
tional, provides high-order derivatives of the DFT energy. 
Note that the if) are arbitrary linear combinations of the 
occupied cigenstates of the Hamiltonian, which are only 
required to minimize the total energy functional. For 
instance, in a crystalline solid t/} could represent Wan- 
nier functions if they are used instead of Bloch functions 
to describe the electronic states. In alternative to the 
path followed here, the 2n + 1 theorem can be demon- 
strated also for constrained functionals, with Lagrange 
multipliers used to impose the orthogonality of the or- 
bitals (Gonze 1995a). 



2. Nonlinear susceptibilities 

Within DFT, the third-order derivatives of the energy, 
Eq. (11), are (Gonze and Vigneron 1989): 

N/2 



E^ = 2j2(^n\AV SCF - Ae n \A?p n ) 

d 2 v 



n=l 

N/2 

+ 2]T ((av„ 

71=1 

N/2 



ox 2 



Ipn ) + C.C. 



+ 2 E(^^^) + ^/^ (3) (n,r 2 ,r 3 ) 

n=l ' ' 

x An(r 1 )An(r 2 )An(r 3 )dr 1 dr 2 dr 3} (122) 

where is the third order functional derivative of the 
exchange and correlation energy with respect to the den- 
sity. In this section, we have dropped the indication of 
the order in A from Aip n which is always the first order 
term. 

The solution of Eq. (30) yields the projection on the 
conduction manifold of the first order change of the wave- 
functions. Therefore we need to recast Eq. (122) in a 
form which does not depend on the projection of \Atpi) 
on the valence manifold. This expression for E^ exists 
since the total energy functional is invariant with respect 
to a unitary transformation within the manifold of the 
occupied orbitals. The required transformation has been 
carried out by Debernardi and Baroni (1994) and by Dal 
Corso and Mauri (1994). After some algebra one obtains: 

N/2 

E^ = 2^<AV>„|P C A1/ SCF P C |AV„> 

n=l 
N/2 

-2 (A^„|P c |A^ ro )(V> m |AV SCF |^„) 



n,m—l 
N/2 



71=1 ^ ^ 



d 2 v 



dx 2 



Ipn ) + CC. 



N/2 



+ 2j2Un 



71=1 



d 3 V 



dX 3 



Ipr. 



i J K^(v u r 2 ,r 3 ) 



xAn(ri)An(r 2 )An(r3)dridr 2 (2r3. (123) 



When the perturbation A is an atomic displacement, as 
in Sec. II.D, E^ gives the first anharmonic corrections 
to the energy. These corrections are responsible, for in- 
stance, for the decay of the phonon modes into vibrations 
of lower frequency. The line width of the phonon lines in 
Raman scattering, after subtraction of isotopic and inho- 
mogeneous broadening, is proportional to E^ if higher- 
order processes are neglected. The comparison of theo- 
retical and experimental values is described in Sec. V.F. 
The generalization of Eq. (123) to metals can be done 
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by the techniques introduced in Sec. II. C. 4. The first 
application has been recently presented (Lazzeri 1999). 

When the perturbation A is an electric field, as in 
Sec. II. C. 2, is proportional to the nonlinear opti- 

cal susceptibilities of a material at low frequency. Un- 
fortunately, in this case, Eq. (123) cannot be directly 
used to compute E^ 3 \ In fact it contains a term 
(V>m|A^V SCJ? |V'n) that becomes ill defined when n = in 
(Eq. (37) cannot be applied: the energy denominator 
vanishes). The same problem arises also when Eq. (123) 
is generalized to mixed third order derivatives of the en- 
ergy and one perturbation is an electric field. These 
derivatives allows one to account, theoretically, for the 
intensities of non-resonant Raman lines (Baroni and 
Resta 1986a) or for nonlinear infrared absorption. 

A solution to this problem has been proposed by Dal 
Corso and Mauri (1994), switching to a Wannier rep- 
resentation for the electronic orbitals, and applying the 
2n + 1 theorem to a total energy functional for a periodic 
insulating solid in a finite electric field originally proposed 
by Nunes and Vanderbilt (1994). This functional, which 
exploits the properties of localized orbitals, is written as: 

JV 

E[{w ,m,},Eo]='^2 E ( w o,m\H + eE x\w^ n ) 

I m,n— 1 

x(2<W,/n - (Wl,n\ w 0,m)) , (124) 

where H is the unperturbed Hamiltonian of the solid, 
Eo is the electric field, x is the position operator, e is 
the electron charge, and {wi, m } are functions — in general 
non-orthonormal — localized around the unit cell iden- 
tified by the Bravais lattice vector, Ri. The function 
wi tm is obtained by translating the function centered 
at the origin by a vector Ri: wi^ m {x) — wo^ m (x — R{). 
In practical applications, wo. m is constrained to van- 
ish outside a localization region of radius R c centered 
at the origin. For simplicity in Eq. (124) we have lim- 
ited ourselves to one-dimensional systems of non interact- 
ing electrons. The electron-electron interaction — when 
dealt with within any self-consistent-field scheme — does 
not yield any additional problems. We stress here that 
the expectation value of x is well defined for any finite 
cut-off radius R c . Furthermore we note that even if no 
orthogonality constraints are imposed on the item's, at 
the minimum they become approximately orthonormal 
as shown by Mauri et al. (1993). 

In analogy with Eq. (123) one can show that the non- 
linear optical susceptibility is given by: 

x (2) = ^ E (V /E l = 



- E ^^o^m^pM 1 ^), ( 125 ) 

m.n—l I ' 

the wW's are solutions of a linear system similar to 
Eq. (46) which can be obtained from the condition 

5E<V /8wW = 0: 

-P c eE x|w ,m) = HP c \w^ m ) 

N 

-EE p >K>Kni#K, m >, (126) 

71=1 I 

where P c = 1 - X)n=i \ w i,n)(wLn\ is the projector on 
the conduction bands in the Wannier representation. 

Eq. (125) is difficult to implement since it requires a 
electronic structure code which calculates localized Wan- 
nier functions. However, Eq. (125) can be rewritten in 
terms of Bloch functions. We recall that the Wannier 
functions are defined in terms of the Bloch functions 
i/)*(x) as: 

£1 f 

W0,n( X ) = 7T dk ^n(x), (127) 

2tt Jbz 

where the integral is done over the first BZ, O is the 
dimension of the unit cell, the Bloch functions are nor- 
malized on the unit cell, and tp k+G (x) = tpn(x); here 
G is a reciprocal lattice vector. Inserting this defini- 
tion in Eq. (125) and using the relationship: xip^(x) — 
~M^n( x ) + e%kx M u n( x )i where u k n {x) = e~ ikx ^(x) are 
the periodic parts of the Bloch wave-functions, one finally 
obtains: 




(128) 



, _fe(i) D fc(i) 
where u n = R c Un ■ 

The formalism, as introduced in this paper, allows 

one to access nonlinear phenomena in a frequency range 

where one can assume that the electrons are in the ground 

state. The generalization to finite frequencies, starting 

from time dependent DFT has been explored by Dal 

Corso et al. (1996). 

III. IMPLEMENTATIONS 

A. Plane waves and pseudo-potentials 

The first — and still today by far the most numerous — 
implementations of DFPT were based on the plane- wave 
(PW) pseudo-potential (PP) method (Pickett 1989). 
PWs have many attractive features: they are simple 
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to use, orthonormal by construction, unbiased by the 
atomic positions. Contrary to what happens with local- 
ized (atomic-like) basis sets, it is very simple to check for 
convergence, by just increasing the size of the basis set, as 
given by the kinetic-energy cutoff. The FFT algorithm 
allows to quickly go back and forth from reciprocal to 
real space. An especially important advantage of PW's 
is the absence of Pulay (1969) terms in the calculation of 
energy derivatives. As a consequence the HF expressions 
for forces and for force constants are valid without any 
correction when a PW basis set is used. 

PW's are used in conjunction with PP's. A PP is a 
fictitious electron-ion interaction potential, acting on va- 
lence electrons only, that mimics the interaction with the 
inner electrons — which are supposed to be frozen in the 
core — as well as the effective repulsion exerted by the 
latter on the former, due to their mutual orthogonality. 
Modern norm-conserving PP's (Hamann, Schliiter and 
Chiang 1979) are determined uniquely from the proper- 
ties of the isolated atom, while the requirement of norm 
conservation ensures an optimal transferability. By the 
latter expression one indicates the ability of the PP's to 
provide results whose quality is to a large extent indepen- 
dent of the local chemical environment of the individual 
atoms. Norm-conserving PP's are angular- momentum 
dependent (i.e. they are non-local operators) and a spe- 
cial care must be taken to ensure that the atomic valence 
(pseudo-) wave- functions associated with them are suffi- 
ciently smooth in the atomic (pseudo-) core, so that they 
can be efficiently dealt with using a PW basis set. Ex- 
perience has shown that the use of PP's is practically 
equivalent to the frozen core approximation within an 
all-electron approach. The pseudo-potential approxima- 
tion implicitly assumes that the energy functional is lin- 
ear with respect to the partition of the total charge into 
core and valence contributions. In some atoms (such as, 
e.g., alkali atoms) the loss of accuracy due to the neglect 
of nonlinearity in the exchange-correlation energy func- 
tional can be intolerably high. For such cases the non- 
linear core correction of Louie et al. (1982) turns out to 
be very useful. 

From a computational point of view, it is very conve- 
nient to recast the angular-momentum dependent part 
of PP's into a sum over a few projectors (Klcinman and 
Bylander 1982). This is called the separable form of PP's. 
The use of PW's and of separable PP's, together with 
FFT and iterative diagonalization or minimization tech- 
niques, allows a fast and efficient solution of the KS equa- 
tions for systems containing up to hundreds of atoms in 
the unit cell. The technical aspects of the implementation 
of the KS equations in a PW-PP framework have been 
extensively described in the literature (see e.g. Pickett 
(1989), Payne et al. (1992), Giannozzi (1995)). 

The implementation of DFPT in a PW-PP framework 
is a straightforward extension of the implementation of 



the KS equations. The only modifications to the theory 
expounded in the preceding chapter are related to the 
non-local character of PP's. In DFT calculations this is 
accounted for by modifying the electron-ion interaction 
term in the energy functional, Eq. (11), as follows : 

E[n(r)\ = F[n(r)} 

N/2 

+2^T / r n {v)V(vy)^ n (v')dvdv', (129) 

where V(r, r') is a sum of atomic nonlocal PP's. The 
results of Sec. II. D must be generalized accordingly. Eqs. 
(21) and (22) become: 



n—1 x 



dVx 



(130) 



and 



8 2 E 
dXidXj 




(131) 



The expression (86) for the electronic contribution to 
force constants is modified as follows: 



el c:f (q) = 




c.c. 



(132) 



The procedure outlined above can be easily extended 
to the case of PP's with non- linear core corrections (Dal 
Corso et al. 1993a). In this case, the exchange-correlation 
functional E xc [n(r)] must be replaced in Eq. (12), by 
E xc [n(r)+n c (r)], where n(r) is the atomic valence charge 
density, and n c (r) is the core charge, or a suitable smooth 
approximation to it (Louie et al. 1982). The exchange- 
correlation potential v xc (r) of Eq. (13) is replaced by 
v xc {n{r) + n c (r)). Energy derivatives will contain addi- 
tional terms. The following contribution must be added 
to first derivatives, Eq. (21): 



dE c 



= J v xc (n(r) + n c (v)) dn ^ dr. 



(133) 



Let us specialize to the case of monochromatic atomic 
perturbations. The following term must be added to the 
screened potential used in the calculation of the linear 
response, Eq. (36): 



0u?(q) 



/ 



dv xc (n) 



dn 



n=n(r)+ri c (r) 



dn c (r) 

0u?(q) 



dr (134) 
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where the core charge n c (r) is written as a sum of ionic 
terms, as in Eq. (87). The force constants, Eq. (86), will 
contain additional terms : 



(q) 



Jv x 



: (n(r) +n c (r)) 



£>V(r) 



au;°( q )a«f( q ; 



■dr 



d(ra(r)+n c (r))\* 3V cc (r) 



0U?(q) 



9 U f(q) 



dr 



(135) 



The matrix elements of the relevant quantities between 
PW's are in the appendix. 



B. Ultra-soft pseudo-potentials 

In typical bulk semiconductors (e.g. Si, Ge, GaAs, 
AlAs) at equilibrium volume 100-150 PW's per atom are 
sufficient for most applications. However, many atoms — 
transition metals, first-row elements like F, O, and to 
a lesser extent N and C — require rather hard PP's to 
ensure transferability, and their treatment demands im- 
practically large PW basis sets. One can try to exploit 
the many degrees of freedom which are present in PP gen- 
eration to obtain softer PP's. Several recipes have been 
proposed to get optimally smooth PP's (for example by 
acting on the form of pseudo-wave-functions in the core 
region). Simple and effective recipes have been described 
by Vanderbilt (1985), Rappe et al. (1990) and Troullicr 
and Martins (1991). 

A more radical approach to the challenge posed by 
hard PP's has been proposed by Vanderbilt (1990) who 
introduced ultras oft pseudo-potentials. In this scheme, 
the orbitals are allowed to be as soft as possible in 
the core regions so that their PW's expansion converges 
rapidly, at the price of giving up the norm conservation 
and the standard ortho-normality constraint of atomic 
orbitals. Orthonormality is recovered by introducing a 
generalized overlap operator which depends on the ionic 
positions. The full electron density is obtained by adding 
to the square modulus of the orbitals an augmentation 
charge localized in the core regions. Despite these tech- 
nical complications, this approach has proved to be ex- 
tremely successful in treating large-scale electronic struc- 
ture problems. 

In the ultrasoft scheme, the atomic PP is separated 
into a local V/ oc and a non-local V NL part. The non-local 
potential is written in the separable form, as a sum of 
projectors. The ionic potential is written as a sum over 
all the atoms I of projectors: 

V NL (r,r>) = £5><° )J tf(r- R,)/?f (r' - R,) (136) 



where the functions (3( (r) and the coefficients D^ 1 are 
computed in an atomic calculation, as described by Van- 
derbilt (1990) and Laasonen et al. (1993). 



The charge density is computed by augmenting the 
square modulus of the orbitals with a term which recovers 
the full valence charge density: 

N/2 

n(r) = 2^>„(r)| 2 

71=1 

N/2 

+ 2 E E E QlA* - R./) Wnl#><# l^> 

n— 1 / ij 
N/2 

= 2^<VgA(r)hU. (137) 

71=1 

Consistently with this definition the orbitals arc sup- 
posed to obey generalized orthogonality constraints 
(tpn\S\4> m } = 5 nm , with an overlap operator S given by: 

S(r,r') = 5{v-r') 

+ EE^(r-R/)A*V-R/) (138) 

I ij 

where q(j = J dr Q\- (r) . The S operator depends on the 
atomic positions and so do the constraints obeyed by the 
orbitals. 

The orbitals are determined by minimizing the total 
energy with the above constraints. This yields a gener- 
alized Kohn-Sham equation: 



H SOF \ip n ) = e n S\ip n ) 
where H„ nF is the KS hamiltonian: 



H, 



2m dr 2 



+ V KS , 



(139) 



(140) 



V KS is the KS potential: V KS = V NL + Vi oc + V Hxc , where 
V Hxc is the Hartree and exchange correlation potential, 
and V NL is the non local part, Eq. (136), in which the 
atomic coefficients D^ 1 are replaced by screened coeffi- 
cients D\- (Laasonen et al. 1993): 



(0)1 



QUv)V eSf {v)dv. 



(141) 



where V e ff = Vi oc + V Hxc . The forces acting on the atoms 
are obtained as in Sec. II. C, differentiating the total en- 
ergy with respect to the atomic displacements and using 
the HF theorem. In the ultrasoft PP case, however, the 
orthogonality constraints change as the atoms move thus 
giving rise to additional terms in the forces. Differentiat- 
ing the generalized orthogonality constraints with respect 
to atomic displacements u™ for the s-th atom we obtain 
the following relationship: 



\S\1>r. 



= -(i>n 



dtp m 

dS 



du° 



V>77 



(142) 
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which, used in the expression of the forces gives: 



N/2 



n=l 



dV K£ 



du° 



dS 



(143) 



where the partial derivative of V KS is: 
dV KS (r 1 ,r 2 ) _ dV NL (r 1 ,r 2 ) 



dul 



dul 



, dVi oc {r 3 ) 

dr 3 — — — A(r 3 ;ri,r 2 ) 



du<i 



+ J dr 3 V eff (r 3 ) 



dA(r 3 ;r 1 ,r 2 ) 



(144) 



with no derivative of the Hartree and exchange and cor- 
relation potential. 

In order to compute the interatomic force constant and 
hence the dynamical matrix we need the mixed second 
derivatives of the total energy with respect to the dis- 
placements Ug and uf of the atoms at sites s and t. 
These expressions have been derived by Dal Corso et al. 
(1997). Here we report only the final results. Taking 
the derivative of the HF forces, one finds that the elec- 
tronic contribution to Cff contains four terms. The first 
one, C^ Q/3 , corresponds to the expectation value of the 
second derivative of the electron-ion potential, and the 
additional second derivative of the overlap matrix: 



C 




d 2 (y NL + v loc ) 



dufdu+ 

d 2 S 
''du?du? 



where the second derivative of V NL is performed at fixed 
charge density as before. The second term C\ 



(2)a/3 



is: 



N/2 r 

= 2 £ 

n=l 



du° 



\P+\4 n )+h.c. 



(146) 



where P c + = 1 — X^m=i S\^m){'4 ) m\ is the projector on the 
conduction-band subspace, h.c. indicates the hcrmitian 
conjugate with respect to the sa, t(3 indices, <p is defined 

as 



\4>L) 



d (V NL + V loc ) ds 



n dul 



Ipn 



(147) 



and again the derivative of V NL is performed at fixed den- 
sity. In the norm-conserving PP scheme, the electronic 



contribution to Cff is simply given by the sum of C. 



(l)a/3 

Hill Ul l_y 

and Cif a0 calculated for S = 1 and A(r) = |r)(r|. In 
the ultrasoft PP scheme one must consider two additional 



contributions to Cff which have no counterparts in the 
norm-conserving scheme. C^"' 3 is the interaction be- 
tween the change of the augmentation charge A t/3 n(r) 
due to the atomic displacement uf (see Eq. (152) below) 
and the change of V Hxc due to the displacement u" (see 
Eq. (27)): 



r (3)af) _ 1 f 



dV Hxc {r) 



A t/3 n(r)dr + h.c. 



(148) 



Finally, C^ a/3 is analogous to C^ a ^ but with the pro- 
jector on the conduction-states subspace replaced by that 
on the valence-state subspace. Since the perturbation 
formalism provides explicitly only P c \^k), the valence- 
state component must be derived from the constraints 
imposed by the orthogonality condition, Eq. (142). One 
finally obtains: 



c (4) a p 



X :> 

' E 

n,m— 1 



OS 



1pm){lpm\ct>l n )+h.C. (149) 



We note that in the norm-conserving scheme, the left 
hand side of Eq. (142) vanishes since S = 1, and Eq. 
(142) allows one to show that the contribution to Cff 
from the valence-states component of \§^t) is zero. In 
the ultrasoft case, Eq. (142) is used to evaluate such a 
component in terms of the unperturbed orbitals, as given 
in Eq. (149). 

The key ingredient to evaluate the dynamical matrix 
>, which can be determined, within first-order 



is P t 



I dip„ 
9u n 



(145) perturbation theory, by solving the linear system: 



(H SCF — e n S)P c 



dlpr, 



where 

~ dV KS 
du? 



dS 



-P+ 



4>r. 



dV K 
du a a 



dS 



Ipr 



(150) 



+ 



/ 



dV H * c (r) 



A(r)|^„)dr. (151) 



Eq. (150) is the generalization to the ultrasoft case of 
the self-consistent linear system given in Eq. (25). The 
perturbing term depends on the variation of the charge 
density <9n(r)/<9w" through dV Hxc (r) / du" . dn(r)/du" is 
a functional of P c \dip n /dUg): 



dn(r) 



N/2 



71 = 1 



^ = 4 E(fel P c + A(r)|^)+A-n(r). (152) 



dii 
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The term A sa n(r), peculiar to the ultrasoft scheme, has 
two contributions: A sa n(r) — S sa n(r) + S sa n ort ho{ r )- 
The former term, 



N/2 

71=1 



0A(r) 



Ipn 



(153) 



accounts for the displacement of the augmentation charge 
at fixed orbitals, whereas the latter, 



8 sa n ortho (v) = 

N/2 



n,m— 1 



OS 



i> m ) <V>m|A(r)|^„}, (154) 
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,r. 



+ $Fjrrfi 



m u m.n 



n.m— 1 



X ( 7/7. 



as 



^ ro )(V ro |A(r)|V„). (155) 
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C. Localized basis sets, all-electron 

All-electron implementations of DFPT based on lo- 
calized basis sets exist for both the Linearized Muffin- 
Tin Orbitals (LMTO) method (Savrasov 1992, Savrasov 
1996) and for the Linearized Augmented Plane Waves 
(LAPW) method (Yu and Krakauer 1994). LMTO and 
LAPW are amongst the most popular all-electron meth- 
ods in DFT calculations. Their extension to DFPT calcu- 
lations is especially useful for systems containing transi- 
tion metals (like e.g. high-temperature super-conductors 
and most ferroelectrics) for which the PW-PP approach 
is not very practical. An earlier implementation using 
a mixed-basis set (localized atomic-like states plus plane 
waves) in a PP formalism (Zein 1992) and a more recent 
one (Heid and Bohnen 1999) are also known. 

Localized-basis-sets and mixed-basis implementations 
are more complex than PW implementations. Part of the 
additional complexity arises from Pulay terms in deriva- 
tives. The origin of such terms is easily understood. The 
first derivative of the energy functional (see Eq. (53)) con- 
tains the HF term, as in Eq. (21), plus a term F, coming 
from implicit dependence through the wave-functions: 



dE f . ^dV(r) , ~ 



(156) 



where 



N/2 



n=l 



3<(r) 
dX 



(H SCF - e n )ip n (r)dr + c.c. (157) 



vanishes if the wave-functions are exact KS orbitals. This 
is not always true if the wave-functions are approximate 
KS orbitals. Let us expand the wave-functions into a 
basis set (f) n , taken to be orthonormal for simplicity: 



Wr)=5>5 n V;(r). 



(158) 



appears because of the orthogonality constraints, similar 
to the Clf 1 " 13 term in the interatomic force constants. 

The generalization of the above formalism to metallic 
systems can be done along the same lines as described 
in Sec. II. C. 4. The presence of the fractional occupa- 
tion numbers modifies the definition of the valence-states 
subspace and the terms S sa n or tho( r ) and C^ Q/3 must be where 
modified accordingly. For instance S sa n ort ho( r ) becomes: 



The solution of the KS equations reduces to the solution 
of a secular equation 



(159) 



H jk = J 4>*{v)H SCF <p k {r)dr. (160) 

By inserting the expansion of the KS orbitals into Eq. 
(157) one finds 

n=l jk ^ 



dc 



-(flj- fc -e n )4 n) +c.c. 



(161) 



The second term vanish exactly (see Eq. (159)). The 
first term does not vanish if the basis set is not complete 
and if the basis set depend explicitly on A. In realistic 
calculations using atomic-centered basis sets, the Pulay 
contribution cannot be neglected. Accurate and reliable 
calculations of forces and of the force constants require a 
very careful account of the Pulay terms, that are instead 
absent if a PW basis set is used. 

LMTO-based linear-response techniques have been 
used to calculate phonon spectra and electron-phonon 
couplings in several elemental metals (Savrasov et al. 
1994, Savrasov and Savrasov 1996) and more recently 
in doped BaBi03 (Meregalli and Savrasov 1998) and 
CaCu02 (Savrasov and Andersen 1996). An extension 
of the method to the calculation of spin fluctuations 
(Savrasov 1998) has been recently done. 

LAPW-based linear-response techniques have been 
first used in the study of the lattice dynamics of CuCl 
(Wang et al. 1994) and later applied to several mate- 
rials: SiC (Wang et al. 1996a), ferroelectrics KNb0 3 
(Yu and Krakauer 1995, Waghmare et al. 1998, Wang 
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et al. 1996b) and SrTi0 3 (LaSota et al. 1998), high-T c 
super-conductor La2Cu04 (Wang et al. 1999). 

The mixed-basis approach has been tested by calcu- 
lating the phonon dispersions for simple metals con- 
taining localized 4d electrons: Ag and Y (Hcid and 
Bohnen 1999), Ru (Heid et al. 2000b), and on sapphire 
(a-Al 2 3 ) (Hcid et al. 2000a). 



IV. OTHER APPROACHES 
A. Dielectric approach 

Historically, the microscopic theory of lattice dynam- 
ics was first formulated in terms of dielectric matrices 
(Pick, Cohen, and Martin 1970). The basic ingredient is 
the inverse dielectric matrix, e _1 (r,r'), relating, in the 
linear regime, the external perturbation AV^ to the to- 
tal electrostatic potential experienced by an external test 
charge: 



AV test (r)= / e- 1 (r,r , )Ay(r , )dr'. 



(162) 



As an alternative the theory can be formulated in terms 
of the electron polarizability, x(r, r'), which gives the 
charge-density linear response to the external perturba- 
tion: 



Ara(r) = J x(r,r')Ay(r')dr'. 



(163) 



These two response functions are simply related as: 

e -1 (r,r')=<S(r-r / )- / -. r x(ri,r)dn. (164) 

J l r_r i| 

Within DFT, one can also define the independent- 
electron polarizability, xo( r , r ')i as the charge-density re- 
sponse to variation of the total KS potential: 

An(r)= J X o(r,r')AV SCF (r')dr'. (165) 

The expression for xo( r > r ') i n terms of KS orbitals has 
the well known form: 



Xo(r,r') = £ 



fn fn 



C(r)V m (r)^(r')^(r') (166) 



where /„ is the occupancy of the state (/„ = 
9((sf —£ n )/ (J ) in the notations of Sec. II. C. 4), and the 
sums over n and m extend to both occupied and empty 
states. As shown in Sec. II. C. 4, only terms involving vir- 
tual transitions from occupied or partially occupied to 
empty or partially empty states contribute. 

Combining Eqs. (163) and (165) and recalling the re- 
lationship between the bare and the KS self-consistent 
perturbing potential, Eqs. (60) and (57): 



Ay SCF (r) = AV(r) + J K(r, r')An(r')dr', (167) 

one gets an integral equation for \'- 
X(r,r') = Xo(r,r') 

+ J Xo{r,r 1 )K(r 1 ,r 2 )x{r2,r')dr 1 dr 2 , (168) 

or equivalently 

X" 1 ^ O - Xo^r, r') - K(r, r'). (169) 

This equation, projected onto a PW basis set, becomes 
a matrix equation, one for each q-point in the BZ, that 
can be solved by matrix inversion. 

The original dielectric matrix approach has the major 
drawback that the perturbation must be described by a 
local potential and thus it cannot be applied to the lattice 
dynamics problem if modern non-local pseudo-potentials 
are employed to describe electron-ion interaction. In fact 
in this case, not only the unperturbed external potential, 
but also the perturbation itself is described by a non local 
operator and Eq. (163) is not appropriate. For these rea- 
sons the calculation of dielectric matrix has been of lim- 
ited utility for the study of vibrational properties — see for 
some early examples the empirical pseudo-potential cal- 
culations by Bertoni et al. (1972), Resta and Baldereschi 
(1981) and Resta (1983) — while it has been successfully 
employed in the study of the macroscopic dielectric prop- 
erties of simple materials (Baroni and Resta 1986b) and, 
more generally, is an essential ingredient in quasi-particle 
GW calculations. In these latter cases, in fact, even if the 
unperturbed external potential is described by non local 
pseudo-potentials, the perturbation of interest is only lo- 
cal and Eqs. (163) and the following ones still apply. 

A modification of the dielectric matrix approach that 
removes its limitation to local pseudo-potentials has been 
devised by Quong and Klein (1992). The response to the 
bare potential is stored into a bare response Arifc(r): 



An fc (r)= / X o(r,r')AV(r')dr' 



(170) 



= E ^^C(r)V m (r)(Vm|Al/|^), 



that can be calculated once for all for non local potentials 
as well. Then Eq. (165) is applied together with Eq. (167) 
to yield An(r): 

An(r) = Anf,(r) 

Xo(r, Ti)K (n, r 2 ) An(r 2 )dndr 2 . (171) 



This equation is solved for An(r) by inverting the kernel 
of the above integral equation (a matrix in a PW basis): 
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(V - r') - J X o(r, ri)A-(n, v')dr^ . (172) 

The second-order change in the energy can then be cal- 
culated as usual. 

This approach has been used for the calculations of 
force constants in Au (Quong 1994), of electron-phonon 
coupling (Liu and Quong 1996) and of thermal expansion 
in metals (Quong and Liu 1997), and of the structural 
stability and electron-phonon coupling in Li (Liu et al. 
1999). 

This modified dielectric matrix approach is conceptu- 
ally rather similar to the original DFPT. The main dif- 
ference with respect to DFPT is the replacement of the 
self-consistency cycle needed in DFPT with the construc- 
tion and inversion of kernel in Eq. (172). This operation 
is rather time-consuming, since it requires the inversion 
of large matrices and an expensive sum over unoccupied 
bands, but must be done only once for any given point in 
the BZ where the vibrations are computed (as opposed 
to the self-consistency in DFPT that must be performed 
for each phonon mode). 

For small systems the overall computational workload 
of the two approaches is similar, but the size of the sys- 
tems that can be treated by the dielectric approach is 
limited by the growing dimension of the kernel operator, 
Eq. (172). 

B. Frozen phonons 

The frequencies of selected phonon modes can be cal- 
culated from the energy differences — or from the forces 
acting on atoms — produced by finite, periodic, displace- 
ments of a few atoms in an otherwise perfect crystal 
at equilibrium. The first such, so called frozen-phonon, 
LDA calculations were performed already in the early 
80's (see for instance Yin and Cohen (1982)). A frozen- 
phonon calculation for lattice vibrations at a generic 
q-vector requires a super-cell having q as a reciprocal- 
lattice vector and whose linear dimensions must be there- 
fore at least of the order ~ 27r/|q|. In practice, the 
size of the super-cell that one can afford to deal with 
has traditionally limited the application of this tech- 
nique to zone-center or selected zone-boundary phonon 
modes in relatively simple materials. However, zone- 
center phonons are also the best characterized because 
they may be Raman- or infrared-active, so that they do 
not require neutron spectroscopy to be detected. More- 
over, the frozen-phonon and linear-response techniques 
may be combined to study anharmonic effects that would 
be otherwise difficult to calculate directly from perturba- 
tion theory (Baroni and Resta 1986a, Debernardi 1999). 

Phonon dispersions along high-symmetry lines in sim- 
ple materials are determined by the so called inter-planar 



force constants (i.e. by the forces acting on planes per- 
pendicular to the phonon wave- vector when another such 
plane is rigidly moved from equilibrium). Lattice vibra- 
tions along some high-symmetry lines in cubic semicon- 
ductors have been calculated in this way using reason- 
ably sized super-cells (Kunc and Martin 1983). More 
recently several authors (Wei and Chou 1992, Wei and 
Chou 1994, Frank et al. 1995, Ackland et al. 1997, Par- 
linski et al. 2000) have started to calculate phonon dis- 
persions from the interatomic force constants determined 
in real space using the frozen-phonon approach. 

The main advantage of the frozen-phonon approach is 
that it does not require any specialized computer code, 
as DFPT does. This technique can in fact be straight- 
forwardly implemented using any standard total-energy 
and force code, and just some care is needed in the eval- 
uation of numerical derivatives. The principal limitation 
is the unfavorable scaling of the computational workload 
with the range of the IFC's, IZifc- In fact, the calcula- 
tion of IFC's within the frozen-phonon approach requires 
the use of super-cells whose linear dimensions must be 
larger than IZifc, thus containing a number of atoms 
N^f ~ 7tf FC . As the computer workload of standard 
DFT calculations scales as the cube of the number of 
atoms in the unit cell, the cost of a complete IFC calcu- 
lation will scale as 3 N at x 1Z^ FC , where N at is the number 
of (inequivalcnt) atoms in the elementary unit cell (the 
factor 3 accounts for the three, generally independent, 
phonon polarizations). The calculation of IFC's using 
DFPT requires instead the evaluation of the dynamical 
matrices on a regular grid of wave-vectors in the Bril- 
louin zone, whose spacing must be chosen of the order 
of the inverse of the range of the IFC's: Aq ~ 2tt/1Zifc 
(see Sec. II. D. 3). The number of q-points in such a grid 
is of the order of TZ^ FC . As the computational cost for 
the calculation of each column of the dynamical matrix is 
of the order of N^ t — and the number of such columns is 
3N at — the total cost for the calculation of the IFC's (and, 
hence, of complete phonon dispersions) using DFPT is of 
the order of Tl\ FC x 3N* t . 

Another problem strictly related to these considera- 
tions is that of the calculation of phonon dispersions in 
polar materials. In Sees. II. D. 2 and II. D. 3 we saw that 
the long-range character of the dipolc-dipole interactions 
in polar insulators determines a non analytic behavior of 
the dynamical matrices as functions of the wave- vector in 
the long wave-length limit. The real-space counterpart 
of this property is that IFC's are long ranged as they 
decay with the inverse cube of the distance. The inter- 
polation of the dynamical matrices in reciprocal space as 
well as the calculation of the long-range tails of the IFC's 
in real space is made difficult by this problem. Within 
DFPT, the standard remedy is the separate treatment 
of the non-analytic part of the dynamical matrix, using 
information on the ionic effective charges and crystal di- 
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electric constant, as explained in Sees. II. D. 2 and II. D. 3. 
Frozen-phonon super-cell calculations, instead, do not di- 
rectly provide such information, 4 that must instead be 
extracted from the limiting behavior of inter-planar force 
constants (Kern ct al. 1999), supplied by a separate cal- 
culation using the Berry's phase approach (King-Smith 
and Vanderbilt 1993, Resta 1994), 5 borrowed from DFPT 
calculations (Parlinski et al. 1998), or fitted to the exper- 
iment (Parlinski ct al. 2000). 

C. Vibrational properties from molecular dynamics 

All the methods described so far are static, zero- 
temperature, methods. In the last 15 years the com- 
bined use of molecular dynamics (MD) and of DFT (Car 
and Parrinello 1985) has become a very powerful tool for 
the ab initio study of condensed-matter systems at fi- 
nite temperature. In MD simulations atomic trajectories 
are generated from the classical equations. Equilibrium 
properties are then estimated as time averages over the 
trajectories which also contain information on the dy- 
namics of the system, i.e. on phonon modes. In fact, 
the vibrational density of states — which exhibits peaks at 
the phonon frequencies — can in principle be computed by 
Fourier transforming the atomic velocity auto-correlation 
function (Rahman 1964). Ab initio MD simulations are 
usually performed using super-cells which contain a small 
number of atoms (from a few tens to a few hundreds) 
and periodic boundary conditions. Because of this, only 
phonons which are zone-center phonons of the super-cell 
are accessible to the simulation. 

The straightforward estimation of phonon frequencies 
from MD simulations suffers from three types of prob- 
lems. First, at low temperature all the systems are 
strongly harmonic and, hence, poorly ergodic. The time 
necessary to reach equilibrium may thus be impractically 
long. In these cases, bringing the system to thermal equi- 
librium may require technical tricks such as coupling to 
Nose-Hoover thermostats (Martyna et al. 1992). 

The second problem is that the simulation time nec- 
essary to attain a frequency resolution, Aw, cannot 
be shorter than r > 2-k/Aui. In practice, this time 
may be too long for first-principles MD simulations. 
This problem can be overcome by using more sophisti- 
cated spectral estimation methods, such as the Multi- 
ple Signal Classification (MUSIC) algorithm (Lawrence 
Marple 1987), to extract information on phonon modes 



4 This is a fact occasionally overlooked in the literature: see 
e.g. Parlinski et al. (1997) and the comment by Detraux et al. 
(1998) 

5 Note that the Berry's phase approach cannot be used to 
calculate the dielectric constant 



from relatively short MD runs. MUSIC exploits the har- 
monic character of phonon modes and has proven to 
be very accurate and useful in simple situations. In 
more complicated cases, where the frequencies are too 
many and too close with respect to the inverse length 
of the run, MUSIC suffers from instabilities that may 
prevent the full determination of the spectrum. An im- 
proved self-consistent MUSIC algorithm has been devel- 
oped (Kohanoff et al. 1992, Kohanoff 1994) to cope with 
instabilities and to extract information on eigen-vectors 
as well. In this algorithm, a first estimation of the fre- 
quencies is followed by the determination of the eigen- 
vectors through a least squares fit of the trajectory in- 
cluding orthogonality constraints. Then, the trajectory 
is projected onto each of the normal modes. At this point 
each projected trajectory contains mainly one frequency 
component. This component is rc-cstimated using MU- 
SIC, and the scheme is iterated until self-consistency in 
the frequencies is achieved. 

Finally when MD simulations arc used to predict the 
temperature dependence of individual vibrational modes 
and the thermal behavior of properties which depend on 
them — the results may depend on the size of the simula- 
tion (super-) cell. In fact, even though in the harmonic 
approximation a mode commensurate with the simula- 
tion cell is strictly decoupled from those which are not, 
this is not the case at high temperature when anharmonic 
effects arc important. As a consequence, the neglect of 
modes which are not commensurate with the simulation 
cell may affect the evaluation of the frequencies of com- 
mensurate modes which, in the harmonic approximation, 
would be directly accessible to the simulation. 

MD simulations are in some sense complementary to 
lattice-dynamical calculations, in the sense that the lat- 
ter are better suited at low temperatures where the for- 
mer are subject to ergodicity problems. Lattice-dynamics 
is by definition limited to the (quasi-) harmonic regime, 
while MD naturally accounts for all the anharmonic ef- 
fects occurring at high temperature, provided the size 
of the simulation cell is large enough to allow a proper 
description of the relevant phonon-phonon interactions. 

V. APPLICATIONS 

A. Phonons in bulk crystals 

Phonon dispersions in crystals have been for a long 
time calculated using model force constants where some 
form of interatomic potentials is assumed and the pa- 
rameters of the model are adjusted so as to reproduce 
some known experimental results. Although this ap- 
proach works reasonably well in some cases it has se- 
rious drawbacks that in many situations call for more 
predictive methods. The typical experimental data used 
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to fix the parameters in the model are the phonon fre- 
quencies themselves. Model force constants can be seen 
as a compact way to encode the available experimental 
input, with very limited predictive power when applied 
to other properties. In many cases only a few selected 
frequencies (usually of infrared (IR) and Raman-active 
modes) are known and the results of model calculations 
in the rest of the BZ are questionable. Even when the 
entire dispersion spectra is known (usually by neutron 
diffraction measurements), the knowledge of phonon fre- 
quencies alone is not sufficient to determine completely 
the force constants: the knowledge of the phonon dis- 
placement patterns would be needed as well. The lat- 
ter experimental information is seldom available and, in 
the few cases in which the pattern were measured, the 
comparison with the results of even the best models is 
sometimes very poor. 

Using DFPT the interatomic force constants are com- 
puted from first-principles and, as they are usually accu- 
rate, both good frequencies and good displacement pat- 
terns are obtained, without need of experimental inputs. 

Most, if not all, calculations described in the following 
sections are done at the LDA level that usually provides 
very good results. A test of the performance of general- 
ized gradient corrections (GGA) in the flavor proposed 
by Pcrdew et al. (1996), has been performed in Favot 
and Dal Corso (1999). It is found that GGA system- 
atically lowers the frequencies of phonon branches with 
positive Grunciscn parameters. This effect is correlated 
with the GGA's expansion of the lattice constant, since 
GGA phonon frequencies computed at the experimental 
lattice constants are higher than the corresponding LDA 
ones. A similar trend is found also for magnetic Fe and 
Ni (Dal Corso and de Gironcoli 2000). In this case GGA 
equilibrium geometries are much superior to LDA ones 
and phonon dispersions are correspondingly closer to the 
experimental results (Dal Corso and de Gironcoli 2000). 
In Diamond, Al, and Cu, LDA and GGA equilibrium ge- 
ometries and phonon dispersions have similar accuracy 
with respect to the experimental data. Si is an excep- 
tion since the LDA phonon dispersions are already in 
very good agreement with experiment and GGA slightly 
worsens the comparison (Favot and Dal Corso 1999). 

In most applications phonon dispersions are computed 
at the theoretical equilibrium geometry (lattice parame- 
ters and internal coordinates). This choice is mandatory 
when the experimental geometry is poorly known, but 
it is also, in our opinion, the most consistent one when 
comparing with experimental data at low temperature. 
Inclusion of thermal expansion may become necessary in 
some cases when comparing with room and higher tem- 
perature data. See V.E for the treatment of thermal 
effects. 



1. Simple semiconductors 

a. Elemental and III-V semiconductors 

The phonon spectra and effective charges of group IV 
semiconductors Si and Ge (diamond structure) and of 
zincblende structure III-V semiconductors GaAs, GaSb, 
AlAs, AlSb were calculated in Giannozzi et al. (1991). 
The calculated phonon dispersions and densities of states 
for the latter four compounds are shown in Fig. 1, to- 
gether with experimental data. Zone-center phonons, 
effective charges, and dielectric constants of nine III-V 
zincblende semiconductors were computed, along with 
their piezoelectric constants, by de Gironcoli et al. 
(1989). The phonon dispersions of Si were calculated 
as well by Savrasov (1996) as a test of the LMTO im- 
plementation of DFPT. Dispersions for InP appear in 
a paper devoted to the (110) surface phonons of InP 
(Fritsch et al. 1995a); dispersions for both GaP and 
InP were published in a study of phonons in GaInP 2 
alloys (Ozolins and Zunger 1998). For all these mate- 
rials, phonon spectra and effective charges are in very 
good agreement with experiments, where available. For 
AlAs — for which experimental data are very scarce — 
these calculations provide the only reliable prediction 
of the entire phonon dispersions. For Si, the calcu- 
lated phonon displacement patterns compare favorably 
to those extracted from inelastic-neutron-scattering ex- 
periments (Kulda et al. 1994). 

In all these materials the interatomic force constants 
turn out to be quite long-ranged along the (110) direc- 
tion. This feature had already been observed in early 
calculations (Kane 1985, Fleszar and Resta 1986) and is 
related to the peculiar topology of diamond and zinc- 
blende lattices, with bond chains propagating along the 
(110) directions. 

The force constants of GaAs and of AlAs are espe- 
cially interesting in view of their use in complex GaAlAs 
systems such as super-lattices, disordered super-lattices, 
and alloys. While the phonon dispersions in GaAs are ex- 
perimentally well characterized, bulk samples of AlAs of 
good quality are not available and little experimental in- 
formation on its vibrational modes have been collected. 
Since several years it has been assumed that the force 
constants of GaAs and those of AlAs are very similar and 
that one can obtain the dynamical properties of AlAs us- 
ing the force constants of GaAs and the masses of AlAs 
(mass approximation) (Meskini and Kunc 1978). The 
DFPT calculations provided convincing evidence that the 
mass approximation holds to a very good extent between 
GaAs and AlAs (Giannozzi et al. 1991). This transfer- 
ability of the force constants makes it possible to calcu- 
late rather easily and accurately the vibrational spectra 
of complex GaAlAs systems (Baroni et al. 1990a, Moli- 
nari et al. 1992, Baroni et al. 1990b, Rossi et al. 1993). 
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Somewhat surprisingly, the mass approximation does not 
seem to be valid when the interatomic force constants for 
a well-known and widely used model, the Bond-Charge 
Model (BCM), are employed. A six-parameter BCM for 
GaAs that gives dispersions that compare favorably with 
experiments and ab initio calculations, yields, when used 
in the mass approximation, AlAs dispersions that are 
quite different from first-principles results. This clearly 
shows that information on the vibrational frequencies 
alone is not sufficient to fully determine the force con- 
stants, even when complete phonon dispersions are ex- 
perimentally available. In order to obtain more realistic 
dispersions for AlAs in the mass approximation, one has 
to fit the BCM for GaAs to both frequencies and at least a 
few selected eigenvectors (Colombo and Giannozzi 1995). 

b. II-VI semiconductors 

The II-VI zincblendc semiconductors ZnSe, ZnTe, 
CdSe, CdTe present some additional difficulties in a PP- 
PW framework with respect to their III-V or group IV 
counterparts. The cation d states are close in energy to 
the s valence states so that the d electrons should be in- 
cluded among the valence electrons. Phonon calculations 
performed several years ago, when the inclusion of local- 
ized d states in the pseudo-potential was difficult, showed 
that the effects of cation d electrons could be accounted 
for also by including the d states in the core and by us- 
ing the nonlinear core correction approximation. The 
results showed an accuracy comparable to that previ- 
ously achieved for elemental and III-V semiconductors 
(Dal Corso et al. 1993a). Similar calculations have been 
more recently performed for hexagonal (wurtzite struc- 
ture) CdS (Debernardi et al. 1997, Zhang et al. 1996) 
and CdSe (Widulle et al. 1999) and compared with the 
results of inelastic neutron scattering experiments. 

c. Diamond and Graphite 

The phonon dispersions of Diamond, together with the 
internal-strain parameter, thermal expansion coefficient 
in the Quasi-Harmonic Approximation, and the mode 
Griincisen parameter dispersion curves, were calculated 
by Pavone et al. (1993). A unique feature found in Dia- 
mond is the presence of an over-bending of the uppermost 
phonon branch, whose frequencies have a minimum at the 
BZ center (instead of a maximum as in the other elemen- 
tal semiconductors Si and Ge). This feature has impor- 
tant consequences for the second-order Raman spectra 
(see paragraph V.A.l.g). 

High-pressure spectra, up to 1000 GPa, for Diamond 
are reported in Xie et al. (1999c). Under pressure, the 
phonon frequencies of the X± and L' 3 modes gradually go 



higher than those of X\ and L' 2 , respectively. The over- 
bending of the uppermost phonon branch decreases with 
the increase of pressure. 

The phonon dispersions of Graphite along the A — T — 
K — M line were calculated in Pavone et al. (1996). The 
dispersions exhibit an over-bending similar to that of Di- 
amond for the in-plane dispersion. The dispersions be- 
tween graphitic planes is very flat. The peculiar behavior 
of low-frequency branches along the T — K line can be 
related to the long extent of IFC's along the zigzag chains 
in the graphitic planes. 

d. Silicon Carbide 

Silicon Carbide (SiC) may crystallize in a large vari- 
ety of tetrahedrally coordinated polymorphs. Phonon 
dispersion curves were calculated for the 3C (cubic 
zincblende), 2H (wurtzite) and 4H hexagonal structures 
(Karch et al. 1994). For the 3C structure, elastic and 
Griincisen constants were calculated as well (Karch et 
al. 1994). The behavior under pressure of phonon disper- 
sions, effective charges, and of the dielectric tensor was 
studied by several authors (Wang et al. 1996a, Wellen- 
hofer et al. 1996, Karch et al. 1996a). The interest 
for the dynamical properties of SiC under pressure was 
prompted by a report (Liu and Vohra 1994) that the 
splitting between LO and transverse optic (TO) modes 
of 6H SiC increases with increasing pressure until P = 
60GPa, then it decreases. This was attributed to a de- 
crease of the effective charges. This interpretation was 
not confirmed by any of the above theoretical studies, 
which pointed out instead that an incorrect volume de- 
pendency for the dielectric tensor was assumed in the 
analysis of experimental data. 

e. Nitrides 

The group-Ill nitrides are very interesting materials 
for optoelectronic applications at short wavelengths, as 
well as in high-frequency and high-temperature elec- 
tronic devices. Several groups performed calculations of 
lattice-dynamical properties for zinc-blende and wurtzite 
GaN (Karch et al. 1998), BN and A1N (Karch and 
Bechstedt 1997), wurtzite A1N, GaN and InN (Bungaro 
et al. 2000), wurtzite GaN (Ruf et al. 2000) and A1N 
(Schwoercr-Bohning ct al. 1999). All calculations yield 
good agreement with available experiments, not only for 
mode frequencies but also for displacement patterns. In 
particular, the calculated eigenvector for the Raman- 
active Ei mode in wurtzite GaN compares well with the 
eigenvector obtained from the study of the isotope ef- 
fect (Zhang ct al. 1997). The comparison with recent 
High- Resolution Inelastic X-Ray scattering measurement 
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in wurtzite GaN (Ruf et al. 2000) and A1N (Schwoerer- 
Bohning ct al. 1999) shows good agreement with scat- 
tering intensities, thus validating the correctness of com- 
puted eigenvectors. 

The phonon dispersions in BN and A1N arc consider- 
ably different from those of III- V semiconductors not con- 
taining first-row elements: phonon dispersions for BN arc 
similar to those of Diamond, while A1N dispersions are 
close to those of SiC. Furthermore, the difference between 
BN and A1N phonon dispersions cannot be explained by 
a simple mass approximation but derives from the quite 
different degree of ionicity and covalent strength of the 
two materials. The marked ionicity of A1N bonding yields 
a pronounced structural and dielectric anisotropy in the 
wurtzite structure, larger than that of wurtzite BN and 
SiC (Karch and Bechstedt 1997). The three-phonon de- 
cay of the LO phonon in two acoustic phonons is not 
allowed in GaN and InN, since the LO frequency is much 
larger than the acoustic frequencies over the entire fre- 
quency spectrum Bungaro et al. (2000). 

The pressure dependence of the dielectric and lattice- 
dynamical properties of both zincblende and wurtzite 
GaN and A1N have been recently calculated in Wagner 
and Bechstedt (2000). 

f. Other semiconductors 

Phonon dispersions of zincblende semiconductor CuCl 
were calculated using the LAPW method (Wang et 
al. 1994). CuCl exhibits large anharmonic effects: in 
particular, many peaks in neutron scattering measure- 
ments disappear at temperatures as low as room temper- 
ature. The calculated phonon dispersions however agree 
well with low temperature experimental results. 

The phonon dispersions for bulk layered semiconduc- 
tor e-GaSe were calculated by Adler et al. (1998). The 
bulk dispersions agree well both with neutron scattering 
results and with surface phonon measurements with in- 
elastic Helium-atom scattering. The calculation of (0001) 
surface phonons at the K point yields small differences 
with respect to the corresponding bulk phonons, as ex- 
pected for a layered material. This rules out previous 
assumptions of anomalous surface phonon dispersions. 
In a subsequent study (Panella et al. 1999) the phonon 
dispersions of the similar material InSe was calculated. 
The bulk dispersions for both GaSe and InSe compare fa- 
vorably with experiments on (001) thin films epitaxially 
grown on hydrogen-terminated Si(lll) (1 x 1) surfaces. 

The face-centered orthorhombic inter-metallic semi- 
conductor AI2RU exhibits strong far-IR absorption by op- 
tical phonons. Zone-center phonon frequencies, effective 
charges and dielectric tensor were calculated in Ogiit and 
Rabe (1996), showing anomalously large Born effective 
charges, in agreement with experiments. The analysis of 
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FIG. 2. Second-order Raman spectrum of AlSb in the Fi rep- 
resentation at room temperature; the solid line represents the 
theoretical calculation, the dashed line is the experimental 
spectrum. (Reproduced from Windl et al. (1996).) 

the valence charge density shows instead that the static 
ionic charges of Al and Ru are negligible. Hybridization 
is proposed as the origin of both of the semiconducting 
gap and the anomalous Born effective charges. 

g. Second-order Raman spectra in simple semiconductors 

Neglecting matrix element effects, second-order Ra- 
man spectra are approximately given by the overtone 
density of states: I(u>) = g(u)/2), where g(u>) is the vibra- 
tional DOS. A more accurate description is obtained by 
combining the ab-initio phonon spectra with phenomcno- 
logical polarizability coefficients. This technique was ap- 
plied to solve the long-standing controversy on the sharp 
peak in the spectrum of Diamond near the two-phonon 
cutoff (Windl ct al. 1993). As already observed in Pavone 
ct al. (1993), the sharp peak is due to a maximum in the 
vibrational density of states, originating from the over- 
bending of the uppermost phonon branch. Such over- 
bending, and the maximum in the vibrational density of 
states, are absent in the other elemental semiconductors 
Si and Ge. Neither two-phonon bound states nor polar- 
izability matrix clement effects are needed to explain the 
peak. 

A similar technique was applied to the calculation of 
second-order Raman spectrum of AlSb in the Ti symme- 
try (Windl et al. 1996). The calculated spectrum, shown 
in Fig. 2, agrees well with recent experimental data. Both 
theory and experiment give clear evidence for the exis- 
tence of an over-bending in the highest frequency branch 
of the phonon dispersion with a saddle point at the T 
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point as predicted by Giannozzi et al. (1991). This evi- 
dence is derived from a general discussion of the critical- 
point behavior. Furthermore, the over-bending in AlSb 
is explained as an effect of the very different masses of 
the Al and Sb atoms in contrast to the over-bending in 
Diamond, whose origin lies in a peculiarity of the force 
constants. 

In the above cases, a one-to-one correspondence be- 
tween the overtone density of states and the second-order 
Raman spectrum is quite visible. This is not the case for 
SiC (Windl et al. 1994). The Ti spectrum of SiC ex- 
hibits three distinct peaks at 1302, 1400, and 1619 cm -1 
which occur in the gaps of the overtone density of states. 
This exemplifies the importance of taking into account 
the coupling matrix elements. 



h. Piezo-electricity in binary semiconductors 

DFPT is used to evaluate the linear response of wave- 
functions to a macroscopic electric field, from which 
using the stress theorem of Nielsen and Martin (1985) — 
one finds the induced stress. Due to the presence of large 
cancellations between contributions of opposite sign, the 
results are very sensitive to the overall accuracy of the 
calculation. Well-converged calculations yield results in 
good agreement with available experimental data in III-V 
(de Gironcoli et al. 1989) and in II- VI compounds, with 
the exception of ZnTe (de Gironcoli et al. 1990). Non- 
linear piezoelectricity in CdTe was studied as well (Dal 
Corso et al. 1993b). It was found that piezoelectricity 
is linear over a wide range of volume-conserving strains, 
while it displays strong nonlincarity whenever the strain 
is not volume conserving. This implies that the observed 
nonlinear effects can be accurately accounted for by the 
linear piezoelectric response of the cubic system at the 
strained volume. 



2. Simple metals and super-conductors 

An early calculation of a phonon dispersions in Nb and 
Mo metals was performed by (Zein 1992), using a mixed 
basis set. The aim of the calculation was to explain the 
presence of a dip in in the (COO) branch of Nb and its 
absence in Mo. 

The accuracy of DFPT for metals (described in Sec. 
II. C. 4) was demonstrated by de Gironcoli (1995) for three 
test cases: face-centered cubic (fee) Al, fee Pb, and body- 
centered cubic (bec) Nb (see Fig. 3). In all cases the 
calculated dispersion curves are in good agreement with 
experiments if an appropriate smearing technique is used. 
Phonons in fee Cu, Ag, Au were calculated as a test case 
for the use of ultrasoft PP's (Dal Corso et al. 1997). The 
alternative linear response technique of Quong and Klein 
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FIG. 3. Calculated phonon dispersions for fee simple metal 
Al and Pb and for the bec transition metal Nb. Solid and 
dashed lines refer to different smearing widths (0.3 and 0.7 eV, 
respectively). Experimental data are denoted by diamonds. 
(Reproduced from de Gironcoli (1995).) 



(1992) was applied to the phonon dispersions and inter- 
atomic force constants of fee Al. Phonons of magnetic bec 
Fe and fee Ni have been recently calculated by Dal Corso 
and de Gironcoli (2000). For these metals, good agree- 
ment with experiments is obtained using a spin-polarized 
GGA functional for the exchange and correlation energy 
(see Fig. 4). Phonons in hexagonal close-packing (hep) 
Ru were calculated and measured with inelastic neutron 
scattering (Heid et al. 2000b). Many phonon anomalies 
have been found. Phonons in fee Ag and in hep Y were 
calculated as a test of mixed basis set technique (Heid 
and Bohncn 1999). 

The main interest of phonon calculations in metals 
is possibly related to transport properties, and notably 
super-conductivity. The availability of accurate phonon 
dispersions in the entire BZ allows the calculation of the 
electron-phonon (Eliashberg) spectral function a 2 F(ui) 
and of the mass enhancement parameter A that enters the 
MacMillan equation for the transition temperature T c to 
super-conductivity. Other important quantities that can 
be calculated are the transport spectral function af r F(uj) 
and the Xtr coefficients, that determine the electrical and 
thermal resistivity in the normal state. In simple metals 
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FIG. 4. Upper panel: calculated GGA phonon dispersions 
(solid lines) for magnetic bcc Fe compared to inelastic neu- 
tron scattering data (solid diamond) and to dispersions calcu- 
lated within local spin density approximation (LSDA) (dot- 
ted lines). Lower panel: calculated GGA phonon dispersions 
(solid lines) for magnetic fee Ni compared to inelastic neutron 
scattering data (solid diamond) and to calculated LSDA dis- 
persions (dotted lines). (Reproduced from Dal Corso and de 
Gironcoli (2000).) 

calculations of A and a 2 F(uj) were performed in Al, Pb, 
Li (Liu and Quong 1996), in Al, Cu, Mo, Nb, Pb, Pd, Ta, 
V (Savrasov and Savrasov 1996), in Al, Au, Na, and Nb 
(Bauer et al. 1998). Transport spectral functions and co- 
efficients (Savrasov and Savrasov 1996, Bauer et al. 1998) 
and phonon line-widths due to the electron-phonon cou- 
pling (Bauer et al. 1998) were also calculated. Fig. 5 
shows the results of Savrasov and Savrasov (1996) for 
a 2 F(uj). 

The calculation of electron-phonon coefficients found 
a remarkable application, beyond simple metals, in the 
study of the behavior of molecular solids S, Se, Te under 
pressure. With increasing pressure these transform first 
to a base-centered orthorhombic (bco) super-conducting 
structure, followed by a rhombohedral /3-Po phase, and 
finally for Se and Te by a bcc phase. 

At the phase transition between the /3-Po and the bcc 
phase, a jump is observed in T c in Te. The origin of 
this jump was clarified (Mauri et al. 1996) through the 
study of phonon dispersions and of the electron-phonon 
interactions. The phonon contribution to the free en- 




Frequency, (THz) Frequency, (THz) 

FIG. 5. (a)-(h) Calculated spectral functions a 2 F(u>) of the 
electron-phonon interaction (full lines) for the eight elemental 
metals considered in Savrasov and Savrasov (1996). The be- 
havior of the electron-phonon prefactor a 2 (uj) [defined simply 
as the ratio a 2 F(lo)/ F(to)] is shown by dashed lines. Symbols 
plots present the results of available tunneling experiments. 
(Reproduced from Savrasov and Savrasov (1996).) 

ergy is shown to be responsible for the difference in the 
structural transition pressure observed in low and room 
temperature experiments. 

In S the /3-Po phase is predicted to be followed by a 
simple-cubic phase that is stable over a wide range of 
pressures (280 to 540 GPa), contrary to what is observed 
in Se and Tc. The calculated phonon spectrum and 
electron-phonon coupling strength (Rudin and Liu 1999) 
for the lower-pressure /3-Po phase is consistent with the 
measured super-conducting transition temperature of 17 
K at 160 GPa. The transition temperature is calculated 
to drop below 10 K upon transformation to the predicted 
simple-cubic phase. 

3. Oxides 

Oxides present a special interest and a special chal- 
lenge for anyone interested in phonon physics. On the 
one hand, many very interesting materials, such as fer- 
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roelectrics and high-T c super-conductors, arc oxides. On 
the other hand, good-quality calculations on oxides are 
usually nontrivial, both for technical and for more fun- 
damental problems. In a straightforward PW-PP frame- 
work, the hard PP of Oxygen makes calculations expen- 
sive: the use of ultrasoft PP's is generally advantageous. 
LDA is known to be not accurate enough in many cases 
(and sometimes the entire band-structure approach is 
questionable in oxides). Many oxides have complex struc- 
tural arrangements. In spite of all these problems, there 
have been several calculations of phonon-related proper- 
ties in oxides from first principles. These calculations are 
described in the remaining of this paragraph, with the ex- 
ception of work on phase transitions that is deferred to 
Sec. V.D. 



a. Insulators 

In alkaline-earth oxides MgO, CaO, and SrO in the 
rocksalt structure, LDA yields good agreement with 
available experimental data for lattice vibrations (Schiitt 
et al. 1994). The investigation of the phonon-induced 
charge-density fluctuations of MgO and CaO at the L 
point of the BZ partially supports the breathing-shell 
model of lattice dynamics and rules out the charge- 
transfer model for this class of materials. Moreover, the 
calculations show that the breathing-like charge-density 
response is more pronounced for oxygen than for the 
cations in these compounds. 

Silicon Dioxide (Si02) is a much-studied prototypi- 
cal tetrahedrally coordinated compound, existing in a 
large variety of different structures. At ambient condi- 
tions the ground state structure of Si02 is a— quartz, 
whose phonon dispersions, dielectric tensor, and effec- 
tive charges were studied by Gonze et al. (1992). Ef- 
fective charges in a— quartz are anisotropic (here calcu- 
lated for the first time in DFPT) and those for O exhibit 
some anomalous character (see following section). Fig. 
6 shows the calculated phonon dispersions for a— quartz. 
IFC's in a— quartz were calculated in a subsequent paper 
(Gonze et al. 1994). The availability of force constants is 
important because it allows to extensively test semiem- 
pirical interatomic potentials that are used in molecular- 
dynamics simulations of silica. Phonons, dielectric ten- 
sor, and effective charges of the high-pressure, sixfold- 
coordinated phase stishovite were studied by Lee and 
Gonze (1994a). 

The phonon dispersions of a— AI2O3 (sapphire) have 
been recently calculated in Heid et al. (2000a), using LDA 
and a mixed basis set. Sapphire has a rhombohcdral 
unit cell containing 2 formula units (10 atoms). A weak 
anisotropy in the dielectric tensor and in the effective 
charges is found. 
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FIG. 6. Phonon band structure of a-quartz along selected 
directions. Symbols: experimental data; lines: theoretical 
results. (Reproduced from Gonze et al. (1994).) 



b. Ferroelectrics 

The phonon frequencies at T, dielectric tensor and ef- 
fective charges of Titanium dioxide (Ti0 2 ) in the rutile 
structure were calculated in Lee and Gonze (1994b). Ru- 
tile Ti0 2 is an incipient ferroelectric: the frequency of 
the TO mode A 2u decreases with temperature but never 
goes to zero. It is found that the Born effective charges 
of Ti0 2 rutile are much larger than the nominal ionic 
charges of Ti (4+) and O (2—) ions (and much larger than 
those of stishovite in spite of the similar structure) . This 
may sound rather counterintuitive but it is typical of all 
ABO3 ferroelectrics in the perovskitc structure (Ghoscz 
et al. 1998b), whose prototypical material is Barium Ti- 
tanate (BaTiOs). The effective charges of BaTi03 and 
of similar compounds exhibit both values that are def- 
initely larger (anomalous)^ than the ionic charge, and a 
strong anisotropy (for O only in the cubic structure: the 
effective charge is anomalous for displacements parallel to 
the Ti-0 bond, normal, that is close to the ionic value, in 
the orthogonal direction). By performing an appropriate 
band-by-band decomposition (Ghosez and Gonze 2000) 
of contributions, this effect can be tracked to the dy- 
namical change of hybridization, mainly between O 2p 
and Ti 3d orbitals (Ghosez et al. 1995). Born effective 
charges for cubic WO3 in the defect-perovskite structure 
(Detraux et al. 1997) and for KNb0 3 (Wang et al. 1996b) 
follow the same pattern. The important role of covalence 
in determining the anomalous polarization was demon- 
strated also by Posternak et al. (1994), by computing 
the effective charges of a fake material, similar to KNb0 3 , 
where covalence was artificially removed by an additional 
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potential. 

Phonons at T for the cubic (ideal perovskite) and for 
the rhombohedral phases of BaTiC>3 were calculated in 
Ghosez et al. (1996). The complete phonon dispersions 
for the cubic structure, together with an analysis of the 
IFC's, can be found in Ghosez et al. (1998a). 

c. High-T c super-conductors 

Although the microscopic mechanism which gives 
rise to super-conductivity in high-T c oxides is still un- 
der active debate, accurate phonon dispersions and 
electron-phonon coefficients constitute an important 
piece of information for understanding the properties 
of these materials. Calculations were performed — at 
the LDA level — for CaCu02 (Savrasov and Andersen 
1996), Bao. 6 Ko.4BiC>3 (Meregalli and Savrasov 1998), and 
La 2 Cu0 4 (Wang et al. 1999). 

In hole-doped (n = 0.24) CaCu02, the phonon dis- 
persions and electron-phonon coupling, for both s— and 
d— wave pairing, were calculated using LMTO linear- 
response techniques. The resulting values of A ~ 0.3 for 
d x 2_ y 2 symmetry and A ~ 0.4 for s symmetry suggest 
that the electron-phonon mechanism alone is insufficient 
to explain the high T c but could enhance another d— wave 
pairing mechanism (Savrasov and Andersen 1996). 

Similar calculations (Meregalli and Savrasov 1998) 
were performed for Bao.6Ko.4Bi03, using the virtual crys- 
tal and mass approximations. The A parameter, in- 
cluding also anharmonic contributions, is found to be 
A = 0.34, a too small value to explain high-T c super- 
conductivity in this system within the standard mecha- 
nism. 

In tetragonal La 2 Cu04, the phonon frequencies and 
eigenvectors were calculated using LAPW linear response 
techniques (Wang et al. 1999). The results (see Fig. 7) 
are generally in good agreement with experiments, with 
the exception of lowest-lying branches involving anhar- 
monic motion of the apical oxygen atoms parallel to the 
CuC>2 planes. The octahedral tilt mode at the X point 
is found to be the most unstable mode throughout the 
BZ, consistent with the observed phase transition to the 
orthorhombic structure at low temperature. The calcu- 
lated dispersion of the highest frequency S3 branch is 
in good agreement with experiment, showing that a pro- 
posed large renormalization of the phonon spectrum by 
a Jahn- Teller electron-phonon interaction is unlikely. 

4. Other materials 

Fullcrcnc Ceo forms a molecular solid that can be doped 
with up to six alkali atoms (K, Rb, Cs) per Cgo- The al- 
kalis lose their electrons that fill the conduction band of 
Ceo, originated by molecular ti a states. Sizable frequency 
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FIG. 7. Calculated phonon dispersion of tetragonal La2Cu04 
along the (£,£,0) (£) and (£,0,0) (A) directions. The fre- 
quencies are in cm -1 and the imaginary frequencies are rep- 
resented as negative numbers. The vertical dashed line cor- 
responds to the boundary of the first BZ. (Reproduced from 
Wang et al. (1999).) 

shifts and a large enhancement in the intensity of the four 
IR-active modes of molecular C 60 is observed when K (or 
Rb, or Cs) is added to solid C 6 o- The ten Raman-active 
modes of molecular C@o also exhibit large shifts. In in- 
sulating bec K 6 C6o the phonons at T and the effective 
charges were calculated (Giannozzi and Andreoni 1996). 
The frequencies of Raman- and IR-active modes are re- 
ported in Table I. It is found that the structural relax- 
ation of the Ceo molecule is primarily responsible for the 
frequency changes, while the change in IR relative in- 
tensities is a consequence of the electron transfer. The 
potassium vibrations are found to lie within the range 
68-125 cm -1 and are well decoupled from the Cgo in- 
tramolecular modes. The DFPT results do not support 
the so called charged phonon model. 

Elemental Boron and Boron-rich solids tends to form 
complex structures formed by assembly of icosahedral 
B12 units. The exact structure and vibrational proper- 
ties of such materials are not well known. The compari- 
son of accurate phonon calculations with IR and Raman 
measurements is of great help in determining the atomic 
structure. 

Icosahedral Boron presents a very sharp peak at 525 
cm -1 whose vibrational character has been for a long 
time denied. New Raman scattering experiments under 
pressure were compared with ab initio lattice dynamics 
calculations (Vast et al. 1997). The very good agreement 
of the mode frequencies and their pressure coefficients 
yields unambiguous assignment of all observed features, 
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TABLE I. K 6 C6o: Theoretical and experimental frequencies 
lj of the optically C6o-like active modes and shifts Aoj with 
respect to Ceo- We use the labeling of the Ih group. The H 9 
modes split into a triplet T 9 (left) and a doublet E 9 (right). 
The experimental values are ascribed according to the calcu- 
lated ordering. Their shifts refer to the centers of gravity. 
Units are cm -1 (Giannozzi and Andreoni 1996). 
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a Data from Eklund et al. (1992) 
b Data from Kuzmany et al. (1994) 
c Data from Martin et al. (1993) 
d Data from Pichler et al. (1994) 



including the 525 cm -1 line which is a highly harmonic 
librational mode of the icosahedron and mainly involves 
bond bending. This mode is also identified in the Raman 
spectrum of other icosahedral boron-rich solids (Vast et 
al. 1997). 

Boron Carbide, B 4 C, is the third hardest material after 
Diamond and cubic BN. The building blocks of the crys- 
tal are distorted icosahedral BuC units, but their pre- 
cise arrangement is still experimentally unknown. The 
structure of icosahedral B4C boron carbide was theo- 
retically determined by comparing existing IR and Ra- 
man spectra with accurate ab initio lattice-dynamical 
calculations, performed for different structural models 
(Lazzari et al. 1999). The examination of the inter- 
and intra-icosahedral contributions to the stiffness shows 
that intraicosahedral bonds are harder than intcricosa- 
hedral ones, contrary to previous conjectures (Lazzari et 
al. 1999). 

Phonons in solid Cl-a typical molecular solid-were 
calculated in Bauer et al. (1995). The calculated in- 
tramolecular distance is too large by 8%, maybe due to 
the pseudo-potentials used, and as a consequence the in- 
ternal phonon frequencies are too low. Intcrmolccular 



distances are in good agreement with experiment 6 and 
so are most of the calculated external phonons (calcu- 
lated at T and Y) (Bauer et al. 1995). 

The phonon dispersions of transition-metal carbide 
NbC were studied as a sample application in a techni- 
cal paper (Savrasov 1996). NbC presents several pecu- 
liar phonon anomalies that are well reproduced by the 
LMTO calculations. 

Zone-center phonons and dielectric properties (too, Z*) 
of cubic rocksalt alkali hydrides LiH and NaH were cal- 
culated in Blat et al. (1991). More complete phonon dis- 
persions for LiH and LiD appeared in Roma et al. (1996). 

B. Phonons in semiconductor alloys and super-lattices 

The calculation of phonon dispersions in systems de- 
scribed by large unit cells or lacking periodicity alto- 
gether presents a special challenge. Disordered systems 
(such as amorphous materials or substitutionally disor- 
dered alloys) can be described in a PW-PP framework 
by periodically repeated fictitious super-cells. However 
the needed computational effort quickly grows with the 
size of the unit cell or super-cell (as oc N a , where N 
is the number of atoms, a ~ 3 4- 4 in practical calcula- 
tions) so that even with the best algorithms and the best 
computers available one is limited to systems having at 
most <~ 100 atoms. Such size may or may not be ade- 
quate for the physical system under investigation. If it 
is not, the brute- force approach must be supplemented 
by a more targeted approach. Typically, accurate first- 
principles calculations in suitably chosen small systems 
are used to set up a computationally manageable model 
for the large system. 

1. GaAs/AIAs super-lattices 

GaAs/AlAs super-lattices and alloys are a very suc- 
cessful example of how first-principles calculations sup- 
plemented by an appropriate model can lead to an accu- 
rate description of the dynamical properties of real sys- 
tems. As already mentioned in Sec. V.A.I. a, the mass 
approximation is the key for obtaining accurate dynami- 
cal matrices for large systems at a modest computational 
cost, once the interatomic force constants in real space 
for GaAs (or for the Ga^A^^As virtual crystal) were 
obtained. 

The interest for phonons in GaAs/AIAs super-lattices 
is a consequence of the progress in epitaxial techniques 



6 The agreement is possibly fortuitous since LDA does not 
treat correctly van der Waals interactions and materials held 
together by them. 
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that allowed the growth of ultrathin super-lattices (with 
a period of < 10 atomic layers), in particular along the 
(100) direction. Owing to the large difference between 
the cationic masses, GaAs and AlAs optic mode occur 
in different frequency ranges. In the super-lattice, op- 
tical vibrations are confined in one or in the other of 
the materials. In perfectly ordered super-lattices, the re- 
lation between confined modes and bulk dispersions of 
the component materials is given by the unfolding model 
(Jusserand and Cardona 1989). 

In ultrathin super-lattice it was not clear however how 
valid the unfolding model was, how much disorder at the 
interfaces (cation intermixing) was present, and what the 
effect of disorder would be. The simulation of the dy- 
namical properties of ordered and of partially disordered 
super-lattices, using the mass approximation and large 
super-cells (Baroni et al. 1990a, Molinari et al. 1992) 
allowed to clarify the problem. It was found that the 
unfolding model is very well verified in ultrathin super- 
lattices, but that a sizable amount of cation intermixing 
is needed in order to explain the details of the Raman 
spectra. These findings later received independent con- 
firmations (Gammon et al. 1991, Kechrakos et al. 1991). 

It should be remarked that the ability to obtain such 
results critically depends on the quality of the modeliza- 
tion used. Empirical models, even good ones like the 
BCM, are too crude and do not adequately describe the 
dispersions of AlAs, and as a consequence, the details of 
the spectrum. 

Calculations of the phonon spectra in more exotic 
GaAs/ AlAs systems proceeds quite in the same way. For 
instance, the vibrational properties of an array of GaAs 
thin wires embedded in AlAs were calculated in Rossi et 
al. (1993). 

2. GaAs/AIAs alloys 

The same approach that was used for GaAs/ AlAs 
super-lattice was also used to clarify the results of Raman 
measurements in GaAs/ AlAs homogeneous alloys. The 
Raman spectra has two distinct peaks, corresponding to 
the vibrations of each cationic species separately. This 
behavior is called two-mode and is typical of all A X B\- X C 
III-V alloys, with the exception of GaaJni-^P. The peaks 
are shifted and asymmetrically broadened with respect to 
the pure materials. This asymmetry was interpreted as- 
suming that Raman-active phonon modes are localized 
on a scale ~ 100 A, but such assumption was challenged 
by other experimental results indicating that phonons 
have well-defined momentum and are coherent over dis- 
tances > 700 A. The results of simulations using a 512- 
atom super-cell (Fig. 8) clearly indicate that the latter 
picture is the correct one (Baroni et al. 1990b). 



3. Si/Ge super- lattices and alloys 

Si/Ge super-lattices and alloys are more difficult not 
only to grow, but also as a subject of theoretical study, 
than GaAs/ AlAs systems. The mass approximation in 
Si/Ge systems is quite poor, yielding errors of as much 
as 20 cm -1 for optic phonons. Moreover in Si/Ge systems 
the lattice parameters of the two components differ by as 
much as 4%, thus giving raise to sizable strain and atomic 
relaxations that must be taken into account. GaAs/ AlAs 
systems instead are almost perfectly matched (their lat- 
tice mismatch is a modest 0.2%). 

In order to achieve the same level of accuracy for Si/Ge 
systems as in GaAs/ AlAs systems, one has to supplement 
the mass approximation with a correction that takes into 
account the effect of strain and of atomic relaxation. This 
goal is achieved by introducing higher-order interatomic 
force constants that are fitted to first-principles results 
for a few selected configurations (de Gironcoli 1992). 
More complex systems can then be simulated by suitable 
super-cells as in GaAs/AIAs. In this way the vibrational 
properties of Si^Gei-^ were studied (de Gironcoli 1992). 
In particular it was shown that this approach is able 
to reproduce the Raman spectra in Sio.sGeo.s, including 
the details due to the local arrangement of atoms, while 
the mean-field approach (Coherent Potential Approxima- 
tion, or CPA) badly fails in reproducing the three-mode 
character of the spectra (de Gironcoli and Baroni 1992). 
The higher-order interatomic force constants were also 
used to study the vibrational properties of ideal and 
realistically intermixed Si/Ge super-lattices (de Giron- 
coli et al. 1993, Schorer et al. 1994, de Gironcoli and 
Molinari 1994). 

4. AIGaN alloys 

The zone-center vibrational properties of wurtzite 
Al x Gai_ x N alloys, over the entire range of composition 
from pure GaN to pure A1N, were studied by Bungaro 
and de Gironcoli (2000) using the mass approximation 
and the arithmetic average of the interatomic force con- 
stants of the two pure materials previously calculated by 
Bungaro et al. (2000). While some of the alloy modes 
display a two-mode like behavior, do not preserve a well- 
defined symmetry and have a large broadening, the LO 
modes, instead, display a one-mode behavior and have 
a well defined symmetry, small broadening, and a pro- 
nounced dependence of the frequency upon alloy compo- 
sition. Therefore these modes are proposed as the best 
candidates for the compositional characterization of the 
alloy (Bungaro and de Gironcoli 2000). 
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FIG. 8. Spectral densities of states of Gao.5Alo.5As along the F — X direction. The position of the peaks are indicated by the 
solid lines in the oj-q plane. (Reproduced from Baroni et al. (1990b).) 



5. GaP/lnP alloys 

A different approach to study phonons in semicon- 
ductor alloys uses suitably chosen small super-cells, the 
special quasi-random structures (Zunger et al. 1990), to 
simulate a disordered system. The evolution of the vibra- 
tional properties in GaP /InP systems with long-range or- 
der was studied using such approach (with a 16-atom cell) 
to calculate the phonon spectra of random Gao.5Ino.5P 
(Ozolins and Zunger 1998). The phonon spectra of pure 
GaP, InP, and of CuPt-type ordered GaInP2 were calcu- 
lated for comparisons. It was found that ordered GaInP2 
and Gao.5Ino.5P have qualitatively different phonon spec- 
tra: ordered GaInP2 exhibits a two-mode behavior, with 
two GaP-like and two InP-like phonon modes, while dis- 
ordered Gao.5Ino.5P exhibits a pseudo-one-mode behav- 
ior: two LO modes, one of GaP and another of mixed 
GaP/lnP character, appear, while the TO modes of GaP 
and InP have merged into a single alloy mode. This is 
in remarkable agreement with experiments (Ozolins and 
Zunger 1998). 

6. Localized vibrations at defects 

Localized vibrational modes of impurities contain a 
wealth of information on the local structure of the de- 
fect. Their analysis require an accurate knowledge of the 
phonon spectra of the host crystal. The interpretation of 
the isotopic fine structure of substitutional impurities in 
III-V semiconductors was performed using DFPT for the 
bulk crystal and a Green's function technique, with re- 
sults that are far superior to those obtained using model 
calculations (Robbie et al. 1995). With these techniques, 



the host isotope fine structure of 12 C:As and 11 B:As lo- 
cal modes in GaAs (Robbie et al. 1995), of the As:P gap 
mode (Grosche et al. 1995) and of B:Ga gap and local 
modes in GaP (Robbie et al. 1996) were successfully an- 
alyzed. 

C. Lattice vibrations at surfaces 

Most DFPT calculations on surface phonons focussed 
on semiconductor surfaces (Fritsch and Schroder 1999, 
Eckl et al. 1997b, Stigler et al. 1998). Calculations were 
performed for GaAs(llO) (Fritsch et al. 1993), InP(llO) 
(Fritsch et al. 1995a), GaP(llO) and InAs(llO) (Eckl et 
al. 1997a), InSb(llO) (Buongiorno Nardelli et al. 1995b); 
for Si(lll) (Ancilotto et al. 1991); for Si(001) (Fritsch 
and Pavone 1995) and Ge (001) (Stigler et al. 1997); 
for H-covered (110) surfaces of GaAs, InP (Fritsch et 
al. 1995b), and of GaP, InAs (Eckl et al. 1997b); for 
H- and As-covered Si(lll) (1 x 1) surfaces (Eckl et 
al. 1997b, Honke et al. 1996b); for Ga- and B-covered 
Si(lll) (V3 x x/3) R30 (Eckl et al. 1997b); for Ge on 
GaAs(llO) (Honke et al. 1996a); for As-covered (110) 
surfaces of GaAs, GaP, InAs, InP (Fritsch et al. 2000). 
The work on semiconductor surface phonons has been de- 
scribed in a recent extensive review article (Fritsch and 
Schroder 1999). In Fig. 9 we show an example of the ac- 
curacy that can be reached by these calculations in III-V 
surfaces. 

Other DFPT calculations were performed for H- 
covered and clean W (110) surfaces (Bungaro et al. 1996), 
for Be(0001) (Lazzeri and de Gironcoli 1998) and for 
Ag(lll) (Xie et al. 1999b) surfaces. The description of 
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FIG. 9. Phonon dispersion of InP, GaP, In As, and 
GaAs(llO). The large shaded area represents the surface pro- 
jected bulk states. Surface-localized and resonant modes are 
indicated by solid and dashed lines. The squares and triangles 
represent data from HREELS and HAS, respectively. Earlier 
LDA-results are indicated by circles. (Reproduced from Eckl 
et al. (1997b).) 



D. Soft phonons and pressure-induced lattice 
transformations 

Many phase transitions, both induced by pressure and 
by temperature, are driven by a lattice instability. This 
may be an elastic instability (leading to a change of shape 
of the unit cell) or a phonon instability (a soft phonon, 
whose energy goes to zero). When soft phonons have a 
nonzero k vector (usually at the border of the BZ) the 
distortion that sets up causes an increase of the dimen- 
sion of the unit cell. The identification of the soft phonon 
responsible for the phase transition is only the first step in 
understanding the phase transition. The next step is usu- 
ally the construction of a realistic model that takes into 
account all the relevant anharmonic interactions respon- 
sible for the stabilization of the low-symmetry structure 
(these may include coupling with the strain, multiple soft 
modes, and so on). 

The earliest calculation using DFPT of a phonon in- 
stability was performed on narrow-gap semiconductors 
GeTe, SnTe, PbTe (Zein et al. 1992). The first two ex- 
hibit a transition for a high-temperature cubic rocksalt 
phase to a low-temperature rhombohcdral phase, respec- 
tively at T ~ 700K and T ~ IAQK. The calculation of 
the dielectric and zone-center phonons in the cubic phase 
yields negative values for u>to m GeTe and SnTe and an 
anomalously low value in PbTe, consistently with exper- 
imental observations. 



the two latter works is deferred to Sec.V.E.3. 

The hydrogen-covered (110) surface of W exhibits 
phonon anomalies, clearly caused by H adsorption, whose 
nature is still unclear. An anomaly in the upper Rayleigh 
phonon branch along the (001) direction is observed 
both in Helium-atom scattering (HAS) and in electron 
energy-loss (EELS) spectra. In the lower branch, a 
similar anomaly is observed by EELS, whereas a much 
deeper one is only detected by HAS. DFPT calculations 
(Bungaro et al. 1996) yield excellent agreement with both 
HAS and EELS for the upper anomalous branch, and 
with EELS for the lower branch, provided that a careful 
sampling of the surface BZ is performed. Such anomalies 
are interpreted as due to Fermi-surface nesting (Kohn 
anomaly). However the calculations do not predict the 
deep anomaly in the lower branch observed by HAS, 
whose nature is still not clear. As a byproduct of these 
calculations, the phonon dispersions for bulk bcc W and 
for the clean (110) surface are calculated, and both are 
found in excellent agreement with experiments (Bungaro 
et al. 1996). 



1. Ferroelectrics 

The ferroelectric transition in perovskite materials, 
whose most famous example is BaTiOa, is closely related 
to a lattice instability. In ferroelectrics (and in general 
in temperature-induced phase transitions) anharmonicity 
plays a fundamental role: harmonic calculations gener- 
ally yield a negative frequency for the soft mode at zero 
temperature. At higher temperature, anharmonicity sta- 
bilizes the soft mode. Accurate phonon calculations are 
in any case the starting point to construct an effective 
Hamiltonian for the ferroelectric transition through the 
use of a localized, symmetrized basis set of lattice Wan- 
nier functions (Rabe and Waghmare 1995). These are 
the phonon analogous of electronic Wannier functions for 
electrons. Furthermore, the mapping of the phonon in- 
stabilities in the full BZ gives a real space picture of the 
ferroelectric instability, even when it involves coordinated 
atomic displacements in several unit cells. 

Ferroelectric phase transitions were studied using 
DFPT in Barium Titanate, BaTi0 3 (Ghosez et al. 1997), 
Potassium Niobate, KNb03 (Yu and Krakauer 1995, 
Wang et al. 1996b), Strontium Titanate, SrTi0 3 (LaSota 
et al. 1997), Lead Titanate, PbTi0 3 (Waghmare et 
al. 1997), Lead Zirconate, PbZr0 3 (Ghosez et al. 1999). 
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In KNbC>3 phonons at the T point, effective charges 
and dielectric tensor were calculated for the cubic, tetrag- 
onal and rhombohedral perovskite structure. In the hy- 
pothetical cubic structure soft modes are present, one of 
these mode is stabilized in the experimentally observed 
tetragonal structure, and all of them are stable in the 
rhombohedral structure that turns out to be the most 
stable one. An effective Hamiltonian is constructed with 
lattice Wannier function by (Waghmare et al. 1998). For 
the cubic structure a complete mapping of the phonon 
dispersion in the Brillouin zone has been computed. The 
results show a soft mode dispersion that exhibits an in- 
stability of a pronounced two dimensional nature in re- 
ciprocal space and suggest a one-dimensional chain type 
instability oriented along the (100) direction of displaced 
Nb atoms(Yu and Krakauer 1995). Similar phonon dis- 
persions and chain type instability were also found in 
BaTi0 3 (Ghosez et al. 1998a). 

SrTiOs exhibits both ferroelectric 

and anti-ferrodistortive instabilities. In the cubic struc- 
ture, phase space of the ferroelectric instability is greatly 
reduced compared to KNbC>3. Anti-ferrodistortive in- 
stabilities exist in one-dimensional cylindrical tubes ex- 
tending along the entire R— M— R line in the BZ. The 
one-dimensional character of these tubes corresponds to 
real space planar instabilities characterized by rotations 
of oxygen octahedra (LaSota et al. 1997). 

In PbTiC>3, phonons at T, R, X, M for the cubic struc- 
ture are used to construct lattice Wannier functions for 
an effective Hamiltonian. In contrast with the results 
for BaTi0 3 and KNb0 3 , a significant involvement of the 
Pb atom in the lattice instability is found. Monte Carlo 
simulations for this Hamiltonian show a first-order cubic- 
tetragonal transition at 660 K. The resulting temperature 
dependence of spontaneous polarization, c/a ratio, and 
unit-cell volume near the transition are in good agree- 
ment with experiment. Both coupling with strain and 
fluctuations are necessary to produce the first-order char- 
acter of this transition (Waghmare et al. 1997). 

The full phonon dispersion relations of PbTiC>3 and 
of PbZrC>3 in the cubic perovskite structure were com- 
puted and compared with previous results for Barium 
Titanate in Ghosez et al. (1999). The comparison (see 
Fig. 10) shows that the change of a single constituent has 
a deep effect on the character and dispersion of unstable 
modes, with significant implications for the nature of the 
phase transitions and the dielectric and piezoelectric re- 
sponses of the compounds. The unstable localized ferro- 
electric mode of PbTiC>3 has a much weaker dispersion 
with respect to BaTiC^. As a consequence the ferroelec- 
tric distortion is almost isotropic in real space. Further- 
more, there is an anti-ferrodistortive instability at the R 
point, not present in BaTiC>3 or KNbC>3. In PbZr0 3 the 
ground state is anti-ferroelectric and is obtained by freez- 
ing mainly modes at R and S. The phonon dispersions 



show therefore even more pronounced instabilities. The 
unstable branches are dominated by Pb and O displace- 
ments. Examination of the interatomic force constants in 
real space for the three structures PbTi03, PbZr03 and 
BaTiC-3 shows that while most are very similar, it is the 
difference in a few key interactions which produces the 
observed changes in the phonon dispersions. This sug- 
gests the possibility of using the transferability of force 
constants to predict the lattice dynamics of perovskite 
solid solutions. 

2. Pressure- induced phase transitions 

The subject of pressure-induced phase transitions has 
become increasingly important since the invention of the 
diamond anvil cell. Beyond a fundamental interest, the 
behavior of matter at very high pressure (such as those 
found in the interiors of the earth and of other planets) 
is relevant in geology and astronomy. 

Contrary to what happens in temperature-induced 
phase transitions, anharmonicity does not necessarily 
play an important role in pressure-induced phase tran- 
sitions. One or more harmonic frequencies may become 
soft as a consequence of the changes in volume and in 
the atomic positions caused by applied pressure. It is 
therefore important to determine if phonon instabilities 
occur and if they occur at lower or higher pressure with 
respect to other possible instabilities. 

a. Cesium Halides 

Cesium Halides — Csl, CsBr, CsCl — crystallize in the 
cubic B2 structure at low pressure. Under high pressure, 
Csl makes a continuous transition to a lower-symmetry 
phase, whose onset is as low as P ~ 15 GPa. The lower- 
symmetry phase was originally thought to be tetrago- 
nal, but it was later identified as an orthorhombic phase, 
approaching an hexagonal closed-packed structure with 
increasing pressure (Mao et al. 1989, Mao et al. 1990). 
CsBr has also a phase transition around P = 53 GPa, 
while it is not clear if such a transition is present in CsCl. 

A detailed study of the phonon spectra and of the clas- 
tic stability of Csl (see Fig. 11) reveals that an elas- 
tic instability leading to a cubic to tetragonal transi- 
tion is in competition with a phonon instability of zone- 
boundary M§ modes (Buongiorno Nardelli et al. 1995a). 
In the framework of Landau theory, a phenomcnologi- 
cal model for the free energy around the transition can 
be set up using the results from first-principles calcu- 
lations. Group theory allows one to restrict the search 
for minimum-energy structures to those whose symme- 
try group is a subgroup of the group of the undistortcd 
structure. The resulting model is still quite complex: 
the order parameter (the amplitude of the soft phonons) 
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FIG. 10. (Color) Calculated phonon dispersion relations of BaTiOs, PbTiOa, and PbZrC>3 along various high symmetry lines 
in the simple cubic BZ. A color has been assigned to each point based on the contribution of each kind of atom to the associated 
dynamical matrix eigenvector (red for the A atom, green for the Batom, and blue for the oxygen). (Reproduced from Ghosez 
et al. (1999).) 
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FIG. 11. Phonon dispersion relation along the £ (110) for Csl 
at equilibrium volume and just above and below the soften- 
ing pressure of the acoustic mode. Negative frequencies 
actually mean imaginary {i.e. negative squared frequencies). 
(Reproduced from Buongiorno Nardelli et al. (1995a).) 



is six-dimensional. The coupling between phonons and 
strain plays a crucial role in favoring in Csl the transition 
to an orthorhombic structure, driven by a soft phonon, 
with respect to the tetragonal one. In CsBr, instead, the 
elastic instability leading to the cubic to tetragonal tran- 
sition may occur before the phonon instability. Finally, 
in CsCl no softness of the zone-boundary phonons and 
a very weak tendency towards the elastic instability is 
observed (Buongiorno Nardelli et al. 1995a). 



b. Cesium Hydride 

Like most alkali hydrides, CsH crystallizes in the rock- 
salt (cubic Bl) structure and undergoes a transition to 
the CsCl (cubic B2) structure under moderate pressure. 
A second transition from the B2 structure to a new or- 
thorhombic phase (assigned to a CrB structure, space 
group -D2D nas been observed in CsH at P ~ 17 GPa 
(Ghandehari et al. 1995). In Saitta et al. (1997) it is 
shown that this first-order transition is intimately related 
to a displacive second-order transition (driven by a soft 
phonon at the M point of the BZ, M^ ) which would oc- 
cur upon application of a shear strain to the (110) planes. 



c. Silicon Dioxide 

Silicon Dioxide exists in a large number of differ- 
ent phases, both crystalline and amorphous. A sur- 
prising and still not fully understood phenomenon (ob- 
served in many other materials as well) is pressure- 
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FIG. 12. Left panel: a-quartz phonon spectra along the T — L 
line for the soft mode at P=30 GPa (solid), P=34 GPa (dot- 
ted), P=43 GPa (dashed). Right panel: dependence of the 
soft-mode frequency at the K point upon pressure. (Repro- 
duced from Baroni and Giannozzi (1998).) 



induced amorphization, taking place for Si02 at room 
temperature with an onset around 15-25 GPa (Kingma 
et al. 1993a, Kingma et al. 1993b). DFPT calcula- 
tions (Fig. 12) show that a zone-boundary phonon, at 
the K point of the BZ, becomes soft at P ~ 30 GPa 
(Baroni and Giannozzi 1998). The relevance of this fact 
to pressure-induced amorphization is still unclear. Ba- 
roni and Giannozzi (1998) suggested that the extreme 
flatness of the acoustic phonon band whose edge goes 
soft may be related to the strange behavior of amor- 
phization. The iso-structural compound AIPO4 bcrlin- 
ite, which undergoes a similar pressure-induced amor- 
phization, has also a soft phonon at the K point of the 
BZ (Baroni and Giannozzi unpublished), and displays 
a similar — although less pronounced — flattening of the 
acoustic band just before amorphization. 

The stability of stishovite (tetragonal rutile structure) 
under high pressure was studied in Lee and Gonze (1997). 
Both a ferro-elastic instability to the orthorhombic CaCl2 
phase and a softening of the B\ g phonon mode at T were 
found. The former precedes the latter, at P = 64 GPa, 
thus leading to a second-order transition to the CaCl2 
structure. 



d. Semiconductors 

In A N B S ~ N octet semiconductors, one would ex- 
pect a sequence of transitions with increasing pressure 
into the structures of more ionic binary compounds: 
zincblende— >NaCl— ■> (3— tin. Earlier experimental and 
theoretical data supported this picture. However newer 
and more accurate angle-dispersive x-ray measurements 
have revealed that the high-pressure structures of most 
III-V semiconductors are more complex than expected on 
the basis of this simple picture. The reason is that the 
NaCl and f3— tin structures may become unstable with 
respect to lattice instabilities (soft phonons). In Ozolins 
and Zungcr (1999) the systematic absence of the NaCl 
phase for all covalent compounds below a critical ion- 
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icity value is explained in terms of the instability of a 
zone-boundary transverse acoustic (TA) phonon at the 
X point, leading to a Cmcm phase. The (3— tin struc- 
ture turns out to be unstable for all but the most cova- 
lent compound due to the presence of a phonon anomaly 
along the c— axis direction in the LO branch (Ozolins 
and Zunger 1999). A similar approach (Kim et al. 1999) 
shows that in most III-V semiconductors the CsCl struc- 
ture, that was believed to follow the (3— tin structure at 
very high pressure, is actually unstable for GaP, GaAs, 
InP, and InAs, due to the softening of TA phonons at 
the M point. A Landau analysis of the phase transition 
leads to two candidate high-pressure structures, respec- 
tively the InBi and the AuCd structure. 

The softening of zone boundary TA modes in Germa- 
nium was studied both experimentally with inelastic neu- 
tron scattering and theoretically with DFPT calculations 
in Klotz et al. (1997). The softening is prevented by an 
unrelated first order transition to the (3— tin structure 
occurring at P = 9.7 GPa. 

e. Solid Hydrogen 

The nature of the high-pressure phases of solid hy- 
drogen is a fascinating subject. In the gas and diluted 
solid phases, quantum effects are very important and 
the homo-polar H 2 molecules, interacting via electric 
quadrupole-quadrupole interactions, freely rotate even at 
low temperature. Upon compression solid Hydrogen un- 
dergoes a first transition to a broken-symmetry phase 
(phase II) where the molecule rotation is frozen. At 
P <~ 150 GPa Hydrogen undergoes a second transition 
from phase II to a phase III, not clearly identified, ex- 
perimentally characterized by a strong increase in IR ac- 
tivity in the vibronic range. A theoretical study combin- 
ing constant-pressure ab initio molecular dynamics and 
density-functional perturbation calculations (Kohanoff et 
al. 1999) addressed this transition and found that phase 
II (identified as a quadrupolar Pca2\ hep phase) becomes 
unstable due to a zone boundary soft phonon. The re- 
sulting candidate structures for phase III have larger unit 
cells, accounting for the many libronic peaks experimen- 
tally observed, and much larger effective charges, leading 
to a strongly increased IR activity, as observed experi- 
mentally. 

f. Metals 

The lattice dynamics of the bec and fee phases of W 
under pressure was studied in Einarsdottcr et al. (1997). 
The bec phase is stable at zero pressure. Under ap- 
plied pressure, the bec phase develops phonon softening 
anomalies for P ~ 1200 GPa. At this pressure, however, 
the fee and hep phases have a lower enthalpy than the 



bec phase. The fee phase of W has elastic instabilities 
at zero pressure that stabilize with increasing pressure 
before its enthalpy becomes lower than that of the bec 
phase. 

Phonon instabilities in bec Sc, Ti, La, Hf were studied 
in Persson et al. (2000). These metals exhibit hexago- 
nal (Sc, Ti, Hf) or double hexagonal (La) closed-packed 
structure at low temperature, while at high temperature 
they become bec. 

3. Other phase transitions 

In the substitutional^ disordered narrow-gap semi- 
conductor Pbi-^Ge^Te, a finite-temperature cubic-to- 
rhombohedral transition appears above a critical con- 
centration x ~ 0.005. A hypothetical ordered cubic 
Pb 3 GeTe4 super-cell is studied as a model for such alloy 
(Cockayne and Rabe 1997). Unstable lattice modes are 
found, dominated by off-centering of the Ge ions coupled 
with displacements of their neighboring Tc ions. A model 
Hamiltonian for this system (using the lattice Wannicr 
function formalism) is constructed and studied via Monte 
Carlo Simulations. A transition temperature of ~ 620 K 
is found for the cubic model, compared to the experimen- 
tal value of <~ 350 K for the alloy. 

E. Thermal properties of crystals and surfaces 

The knowledge of the entire phonon spectrum granted 
by DFPT enables the calculation of several important 
thcrmodynamical quantities and of the relative stability 
of different phases as functions of temperature. The first 
calculation of a thermal property (expansion coefficients 
in Si) using DFPT dates back to 1989 (Fleszar and Gonze 
1989). 

The thcrmodynamical properties of a system are deter- 
mined by the appropriate thermodynamical potential rel- 
evant to the given ensemble. In the ensemble where the 
sample volume and temperature are independent vari- 
ables, the relevant potential is the Helmholtz free energy, 
F = E — TS. For a solid in the adiabatic approximation 
the free energy can be written as the sum of an elec- 
tronic and a vibrational term. The electronic entropy 
contribution is easily evaluated in metals, although usu- 
ally neglected, whereas it is totally negligible for insula- 
tors: F e i ~ E e i. The key quantity to calculate in order 
to have access to the thermal properties and to the phase 
stability is the vibrational free energy F p h- 

Far from the melting point, the vibrational free energy, 
Fph, can be conveniently calculated within the quasi- 
harmonic approximation (QHA). This consists in calcu- 
lating Fph in the harmonic approximation, retaining only 
the implicit volume dependence through the frequencies: 
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FIG. 13. Lower panel: Zero-pressure free-energy (solid lines) 
and internal-energy (dashed lines) curves for the a and (3 
phases of Tin as functions of temperature. The thin verti- 
cal dotted line indicates the theoretical transition tempera- 
ture, while the experimental value for T c is shown by the 
arrow. Ao = 359 cal/mole is the T — K free-energy 
difference — including the zero-point contribution — while A = 
482 cal/mole indicates and the latent heat adsorbed in the 
a <-> (3 transition. Finally, the inset displays the temperature 
dependence of the vibrational entropies of the two phases, 
(reproduced from Pavone et al. (1998).) 

F ph {T, V) = -fcsTlog (Tr e -KMV)/k B T^ ) (m) 

where T~i p h(V) is the phonon hamiltonian at a given vol- 
ume. In terms of the phonon spectra, F p h can be written 
as 

F p h(T, V) = -fc B T]T log (j2^ (n+1/2)h ^ k(V)/kBT ) ■ 

(174) 

Once the sum over occupation numbers n is performed, 
one gets the final formula: 

F ph (T,V) = k B Tj2^g[2smh(hcj lk (V)/2k B T)}. (175) 
i,k 

In practical calculations, the force constants are calcu- 
lated at a few volumes and interpolated in between to 
get the volume dependence. Once the phonon spectra on 
the entire BZ is available, the calculation of F p h reduces 
to a straightforward integration over the BZ. 

The QHA accounts only partially for the effects of the 
anharmonicity, through the volume dependence of the 



phonon spectra (this is clearly an anharmonic effect: the 
perfect harmonic crystal would have no volume expan- 
sion with the temperature) but it turns out to be a very 
good approximation at temperatures not too close to the 
melting point. 

The quantities that can be calculated in the QHA 
include equilibrium lattice parameters and elastic con- 
stants, specific heat, thermal expansion coefficients, as a 
function of the temperature. Corrections due to quan- 
tum fluctuations [zero-point motion) at zero tempera- 
ture can be estimated as well. The comparison of the 
free energies of different phases (not related by symme- 
try relations, unlike those considered in V.D) yields the 
relative stability as a function of the pressure and of the 
temperature. Most calculations in this field have been 
performed on simple systems, but there are a few ex- 
amples of applications to surfaces, notably to anomalous 
thermal expansion. 

1. Metals 

The equilibrium lattice parameter and thermal expan- 
sion coefficients for bcc Li and fee Al and Na were com- 
puted in Quong and Liu (1997). The relative stability 
of the various polymorphs of Li (bcc, fee, hep, 9R) was 
examined in Liu et al. (1999). It is found that the trans- 
formation from the 9R structure at low temperatures to 
the bcc phase upon heating is driven by the large vibra- 
tional entropy associated with low-energy phonon modes 
in bcc Li. The strength of the electron-phonon inter- 
action in Li is calculated and found to be significantly 
reduced in the low-temperature 9R phases as compared 
to the bcc phase, consistently with the observed lack of 
a super-conducting transition in Li (Liu et al. 1999). 

The thermal properties of fee Ag, plus the Griineisen 
parameters, were calculated in Xie et al. (1999a). 

The a— (3 phase transition in Tin was studied in Pavone 
et al. (1998). At T = K the free energy of the (3 
phase lies ~ 359 cal/mole above that of the a structure. 
The narrower frequency range spanned by the vibrational 
band in the (3 phase makes its entropy larger at high tem- 
perature. As a consequence, the free energies of the two 
phases, shown in Fig. 13, equal each other at a temper- 
ature of <~ 38 C, in close agreement with the observed 
transition temperature T c = 13 C (Pavone et al. 1998). 

2. Semiconductors and Insulators 

Calculations were reported of the thermal properties, 
plus the Griineisen parameters, of Si (Rignanese et al. 
1996); of the thermal expansion coefficients of Diamond 
(Pavone et al. 1993); of the thermal expansion coefficient 
and specific heat of cubic SiC in the 3C structure (Karch 
et al. 1994). 
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In Diamond, thermal properties were calculated at 
high pressure, up to 1000 GPa (Xie et al. 1999c). The 
P-V-T equation of states has been calculated from the 
Helmholtz free energy. The thermal expansion coefficient 
is found to decreases with the increase of pressure, and at 
ultrahigh pressure (700 GPa), Diamond exhibits a nega- 
tive thermal expansion coefficient at low temperatures. 

The temperature dependence of the diamond-/3-Tin 
phase transition in Si and Ge, that occurs under pressure 
at ~10 GPa, was calculated in Gaal-Nagy et al. (1999) 
using the QHA. 

Thermal properties (specific heat, entropy, phonon 
contribution to the free energy, atomic temperature fac- 
tors) for Si02 in the a— quartz phase and in the high- 
pressure stishovite phase were calculated in Lee and 
Gonze (1995). 

The ability of the QHA to yield thcrmodynamical 
properties of materials over a considerable pressure- 
temperature regime was asserted in recent papers (Karki 
et al. 1999, Karki et al. 2000). In these works the thcr- 
moelastic properties of MgO was calculated over a wide 
range of pressure and temperature. Thermodynamical 
potentials and several derived quantities (such as the 
temperature dependence of elastic constants at high pres- 
sures) were computed and successfully compared with 
experimental data (Karki et al. 1999, Karki et al. 2000). 

3. Surfaces 

The thermal expansion of surfaces was calculated for 
Ag(lll) (Xic et al. 1999b) and for Be(0001) (Lazzeri and 
de Gironcoli 1998). 

In Ag(lll) the top-layer relaxation changes from an 
inward contraction (-0.8%) to an outward expansion 
(+6.3%) as the temperature increases from T = K to 
1150 K, in agreement with experimental findings. The 
calculated surface phonon dispersion curves at room tem- 
perature are in good agreement with Helium-scattering 
measurements (Xie et al. 1999b). 

At the Be(0001) surface, an anomalously large thermal 
expansion was recently observed in low-energy electron 
diffraction experiments. The calculations were tested in 
bulk Be, where they describe very well the thermal ex- 
pansion, and checked against first-principles molecular 
dynamics simulations for the surface. The large ther- 
mal expansion is not found. The discrepancy could be 
explained assuming that the actual surface is less ideal 
than assumed (Lazzeri and dc Gironcoli 1998). 

4. Alloys 

The vibrational contribution to the free energy is 
known to affect the phase stability of alloys. The impor- 
tance of such effect was examined in two metallic alloys, 



15 

r- 10 
E 

fei 5 




Carbon 


Silicon / 


Germanium 


















■ ■■ 



0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5 
kT/Hm_. 



FIG. 14. temperature dependence of the full width at half 
maximum, 2r, of the zone-center optical phonon in Diamond, 
Si, and Ge. Solid lines are the result of DFPT calculation; 
squares represents experimental data. (Reproduced from De- 
bernardi et al. (1995).) 

the Cu-Au system (Ozolins et al. 1998) and the Re-W 
systems (Persson et al. 1999). 

Cu-Au systems were studied using a combination of 
DFPT and of cluster expansion methods. The vibra- 
tional free energy of the alloy is calculated by a clus- 
ter expansion over a small set of representative ordered 
structures having small super-cells, quite in the same way 
as the configurational free energy is calculated. Anhar- 
monic effects are taken into account through the QHA. 
The results indicate that the vibrational free energy con- 
tributes significantly to the phase stabilities and thermo- 
dynamic functions of the CuAu system. In particular 
it tends to stabilize the compounds and alloys with re- 
spect to the phase-separated state, and lowers the order- 
disorder transition temperatures. It is found out that 
the vibrational free energy, not the vibrational entropy, 
is the relevant quantity, due to the larger thermal expan- 
sion coefficients of the alloy with respect to the ordered 
ground states (Ozolins et al. 1998). 

A somewhat similar approach was applied to the study 
of the dynamical and thermodynamical stability of the 
bec and fee disordered Re^Wi-a system. As a byprod- 
uct, the phonon dispersion curves for fee and bec Re were 
calculated. Fee Re is dynamically stable but has pro- 
nounced phonon anomalies; bec Re exhibits phonon in- 
stabilities in large parts of the BZ, similar to those found 
in fee W. Due to the instabilities in bec Rc and fee W 
the vibrational entropy, and therefore the free energy, is 
undefined. The problem is circumvented by using the 
virtual crystal approximation to calculate the phonon 
dispersions of disordered Re^Wi-^ and by applying a 
concentration-dependent nonlinear interpolation to the 
force constants. A region where the bec phase would be- 
come thcrmodynamically unstable towards a phase de- 
composition into disordered bec and fee phases is found 
(Persson et al. 1999). 

F. Anharmonic effects 

Third-order energy derivatives can be calculated di- 
rectly from the linear response using the 2n + 1 theorem 
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(See Sec. II. F). The first calculation of anharmonic force 
constants using DFPT and the 2n + 1 theorem was per- 
formed for the lifetimes of zone-center optical phonons in 
Diamond, Si, and Ge (Debernardi et al. 1995) due to its 
anharmonic decay into two phonons of lower frequency. 
This is the main contribution to the line-width of such 
Raman-active modes, once isotopic and other inhomo- 
geneous broadening contributions are subtracted. The 
temperature and pressure dependences of anharmonic 
lifetimes were calculated as well (see Fig. 14 for the for- 
mer). The results are in good agreement with experi- 
mental data. The microscopic mechanisms responsible 
for the decay are found to be different for different ma- 
terials and to depend sensitively on the applied pressure 
(Debernardi et al. 1995). 

Subsequent work was performed in zinc-blende semi- 
conductors GaAs, AlAs, GaP, InP (Debernardi 1998) and 
in SiC (Debernardi et al. 1999). In III-V semiconductors 
the line-width of Raman-active modes, both transverse 
and longitudinal, and their temperature dependence were 
computed. For longitudinal phonons a simple approxi- 
mation consisting in neglecting the effect of the macro- 
scopic electric field in the anharmonic terms yields good 
results (Debernardi 1998). In 3C (cubic zincblende) SiC 
the pressure dependence, up to 35 GPa, of the line- widths 
of the LO and TO modes at the BZ center was calcu- 
lated. An anomalous behavior is found: the line-width 
of the transverse mode changes very little with pressure, 
while the longitudinal mode shows a monotonic increase 
up to 26 GPa, decreasing abruptly above this pressure. 
The results are in good agreement with new experimental 
data, up to 15 GPa (Debernardi et al. 1999). 

A different approach to anharmonicity combines 
DFPT to calculate harmonic force constants and 
frozen phonons to calculate higher-order force constants 
through numerical differentiation. This approach is less 
elegant and more sensitive to numerical accuracy than 
the use of the 2n + 1 theorem. However it allows to cal- 
culate quartic anharmonic terms and those third-order 
energy derivatives whose calculation with the 2n + 1 the- 
orem is hindered by technical difficulties, such as Raman 
cross sections in the nonresonant limit (Giannozzi and 
Baroni 1994). Such approach was first used to calcu- 
late the contribution of quantum fluctuations (zero-point 
motion) to the bulk modulus of Diamond, to the dielec- 
tric constant of Diamond, Si, Ge, and to calculate the 
temperature dependence of these quantities (Karch et 
al. 1996b). More recently, the anharmonic shift of the 
Raman frequency of Diamond and Si, in which both cu- 
bic and quartic anharmonic terms are equally important, 
have been calculated (Debernardi 1999). The tempera- 
ture dependence of the Raman frequency and the contri- 
bution of zero-point motion are calculated as well as the 
Raman line-width. The same quantities were calculated 
in Ge using the 2n+ 1 theorem for the cubic terms and as- 



suming that the quartic anharmonic force constants may 
be approximated by those of silicon (Lang et al. 1999). 

G. Isotopic broadening of Raman lines 

In the harmonic approximation, the effect of isotopic 
disorder or of isotopic substitution is only reflected in the 
change of masses in the dynamical matrix. The frequency 
shift caused by isotopic substitution depends on phonon 
displacements pattern and may be used as a probe of 
the latter. DFPT results were applied to the interpre- 
tation of Raman experiments in isotopically substituted 
C 60 (Guha et al. 1994), in GaN(Zhang et al. 1997), and 
in SiC (Widulle et al. 1999). 

Isotopic disorder is also responsible for a broadening 
and further shift of the Raman lines, beyond the virtual- 
crystal approximation which replaces the masses of each 
individual atom with their composition-averaged value. 
Disorder effects are usually treated within some kind of 
mean-field approximation (such as the Coherent Poten- 
tial Approximation, or CPA). In Stcininger et al. (1999) 
both the CPA and a super-cell approach were used to 
calculate the broadening of Raman peaks in Si, Ge, and 
a-Tin, using IFC's calculated from first principles. 

In a recent work (Vast and Baroni 2000) a new method 
to study the effects of isotopic composition on the Raman 
spectra of crystals beyond mean-field approximation was 
presented. The Raman cross section is expressed as a di- 
agonal element of the vibrational Green's function, which 
is accurately and efficiently calculated using the recursion 
technique. The method was applied to Diamond — where 
isotopic effects dominate over the anharmonic ones — as 
well as to Ge, where anharmonic effects are larger. In 
both case the results are in very good agreement with 
experiments (Vast and Baroni 2000). 

Other effects of isotopic substitution appear beyond 
harmonic approximation. The molecular volume of crys- 
tals depends on isotopic masses through the zero-point 
motion. This is a tiny but measurable effect that can 
be calculated using the quasi-harmonic approximation 
(Sec. V.E). The dependence of the crystal lattice con- 
stant on isotopic composition was calculated for elemen- 
tal semiconductors Diamond, Si, and Ge (Pavone and 
Baroni 1994) and for compound zincblende semiconduc- 
tors GaAs and ZnSe (Debernardi and Cardona 1996). 
In the latter case, the temperature dependence of the 
derivatives of the lattice constant with respect to both 
the anion and the cation mass was computed, together 
with the linear thermal expansion coefficients and mode 
Grimcisen parameters. 

The dependence of the Raman line- width in 3C-SiC on 
both isotopic composition and temperature for longitu- 
dinal and transverse modes is calculated in Debernardi 
and Cardona (1998). The line-widths exhibit a marked 
and nontrivial dependence on isotopic composition. 
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H. Vibrational broadening of electronic core levels 

Phonons cause a temperature-dependent broadening 
of core-level spectra. The calculation of such broadening 
requires an accurate description of both phonon spec- 
tra and of core-excited electronic states. A theoretical 
framework to deal with such problem has been provided 
by Mader and Baroni (1997), together with an appli- 
cation to the Is core exciton of Diamond, for which a 
strong vibronic coupling is present. The sudden approx- 
imation (optical absorption and emission taking place in 
much shorter times than typical phonon response times) 
is assumed, leading to a Franck-Condon picture. The 
electronic and vibrational degrees of freedom are consis- 
tently treated using DFT. A suitable pseudo-potential, 
generated on the excited atomic configuration, allows to 
simulate the excited electronic states. Due to the large 
lattice distortion, anharmonicity cannot be neglected. In 
the limit of infinite core-hole lifetime, anharmonic effects 
are included within a self-consistent phonon approach. 
The Stokes shift for the Is exciton level is found to be 
about 3 eV (the harmonic approximation would overes- 
timate it by about 1.4 cV) and the phonon broadening is 
about 2 cV (Mader and Baroni 1997) . 

VI. CONCLUSIONS AND PERSPECTIVES 

The results reviewed in this paper witness the blossom- 
ing of ab-initio lattice-dynamical calculations in solids, 
based on density-functional (perturbation) theory. Our 
ability to predict from first-principles phonon-related 
properties of materials depends both on the accuracy of 
the ab-initio calculation of lattice vibrations, and on the 
quality of the approximations that arc needed to relate 
these calculations to the specific property one is inter- 
ested in (such as, e.g., the electric conductivity or the 
temperature dependence of the crystal volume). The ac- 
curacy of the calculations can be appraised by compar- 
ing the calculated frequencies with infrared, Raman, or 
neutron-diffraction experiments. It is fair to state that 
lattice dynamics is the one field of solid-state physics 
where the accuracy of ab-initio calculations can compete 
with absorption or diffraction spectroscopies. Of course, 
it would be vain to put much effort in the calculation 
of quantities that can be measured with a comparable 
or better accuracy. The real value of first-principles cal- 
culations stems from their ability to provide unbiased 
predictions for those materials and in those cases which 
are not easily accessible to the experiment. 

Even though necessarily over-simplified, the physical 
conditions of the sample being studied numerically are 
under our total control and can therefore be varied at 
will. This allows to assess to quality and validity of 
models that relate the atomistic and electronic structure 



of materials, which is usually unknown, to their macro- 
scopic, and experimentally accessible, properties. Once 
the accuracy of the calculated phonon frequencies has 
been assessed, the agreement of the predictions for de- 
rived quantities gives indications on the validity of the 
approximations used to derive them. In the case of the 
dependence of inter-atomic distances upon temperature, 
for instance, their comparison with experiment gives in- 
dications on the validity of the quasi-harmonic approx- 
imation used to calculate them. Rather surprisingly, it 
turns out that in a variety of cases, this approximation 
gives precise results up to near the melting temperature. 

We conclude by arguing that the field of lattice- 
dynamical calculations based on density-functional per- 
turbation theory is ripe enough to allow a systematic 
application of DFPT to systems and materials of increas- 
ing complexity. The availability to a larger community of 
scientists of the software necessary to perform such cal- 
culations will make them a routine ingredient of current 
research, much in the same way as it has occurred for 
standard DFT calculations over the past five years or so. 
Among the most promising fields of application, we men- 
tion the characterization of materials through the predic- 
tion of the relation existing between their atomistic struc- 
ture and experimentally detectable spectroscopic proper- 
ties; the study of the structural (in)stability of materials 
at extreme pressure conditions; and the prediction of the 
thermal dependence of different materials properties us- 
ing the quasi-harmonic approximation. 
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APPENDIX A: PLANE-WAVE PSEUDO-POTENTIAL 
IMPLEMENTATION 

1. Pseudo-potentials 

Norm- Conserving PP's are relatively smooth func- 
tions, whose long-range tail goes like —Z v e 2 /r where Z v 



46 



is the number of valence electrons. There is a different 
PP for each atomic angular momentum I: 



V = V loc {r) + Y,Vi{r)P u 



(Al) 



where Pi = J2m=-i \l m )(l m \ 1S the projection operator 
on states of angular momentum I, \lm) being the angu- 
lar state with (l,m) angular quantum numbers. Usually 
Vi oc (r) ~ —Z v e 2 /r for large r so that the Vi(r) are short- 
ranged. PP's in this form are called semi-local: 

V(r,r')=V loc (r)6(r-r') (A2) 
+ J2Yi m (r)V l (r)6(r-r')Y* n (r'), 

Ira 

where the Yi m (r) = (r\lm) are the usual spherical har- 
monics. 

PP's are usually recast into the Klcinman and Bylan- 
dcr (1982) separable form. Each PP is projected on to 
the atomic reference pseudo- wave-functions </>f s : 



v = v loc + v L + J2 



r <ciw> 



(A3) 



where Vz,(r) is an arbitrary function, Vj'(r) = Vi(r) — 
Vl(t). By construction the original PP and the projected 
PP V have the same eigenvalues and eigenvectors on the 
reference states (f>f s . The separable form may badly fail 
in some cases due to the appearance of spurious ghost 
states. Some recipes have been devised to avoid such 
problem (Gonze et al. 1990, Gonzc ct al. 1991). 

2. Matrix elements 

A PW basis set is defined as 

(r|k + G) = i e 4 ( k+G )- r , t-\k + G\ 2 <E cuU (A4) 

where the G's are reciprocal lattice vectors, k is the wave- 
vector in the BZ, V is the crystal volume, and E cut is the 
cutoff on the kinetic energy of the PW. 

In the PW-PP implementation, Sec. III. A, we assume 
that the electron-ion potential Vi on (r, r') is written as 

V ion (r, r') = ^ w s (r - R ; - t s , r' - R ; - r s ) (A5) 

is 

where v s is the PP for the s-th atomic species, whose 
general form is 

v s (r,r') = v Si i oc (r)S{r - r') +J2 v sA r ^')- ( A6 ) 

i 

The plane-wave matrix elements of the above operator 
are given by: 



(k + G|« a |k + G / )=«a,ioc(G-G / ) (A7) 
+ 5^t; a ,,(k + G,k + G'), 



where: 

v s j oc (G) = ^ J v St i oc {r)e 

and 

^,z(ki,k 2 ) = - 



iGr 



dr, 



(A8) 



e- tki r w s j(r,r')e ik2 - r drdr' . (A9) 



For PP's in the semi-local form, Eq. (A2): 

v Stl (r,r') = P;/ S ,;(r), (A10) 

and 

^(ki,k 2 ) = ^(21 + l)P ; (ki • k 2 ) (All) 



/>oo 

x / r 2 ji(kir)ji(k 2 r)f Si i(r)dr, 
Jo 



where the j/'s are spherical Bessel functions, the Pi's are 
Legendre polynomials of degree 

For PP's in the separable form, Eq. (A3): 



and 



v.^rj) = c s Mr)(3 Stl (r'), 



u S; ;(ki,k 2 ) = c S; ;/3*,(ki)/3 Si /(k 2 ) 



(A12) 



(A13) 



where the /3(k)'s are the Fourier transform of /3(r), as in 
Eq. (A8). 

The matrix elements between PW's of the derivatives 
of the ionic potential, Eq. (88), are 



k + q+G 



dV lc 



k + G') = 



3u?(q) 

~ i ( qa + G a -G> a )e-^ +G - G '> T ° 
x (V s ,; oc (q + G- G') 

+ ^^(k + q+G,k + G')). 



(A14) 



The screening contribution to dV S cF/du"(q), which is a 
local potential in DFT, can be advantageously evaluated 
in real space and transformed back to reciprocal space by 
FFT techniques. 

The matrix elements of the commutator [H SCF ,r] be- 
tween PW's are: 



(k 1 |[H SOP ,r]|k 2 > = — k 1( 5 klk2 
m 



(A15) 



d d 
9ki <9k 2 



_i^; e -*(ki-k a )-T. (__ + __) ^ (kl ,k 2 ), 
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where 5kik 2 = 1 if ki = k2, <5kik 2 = otherwise. 

The matrix elements of the second derivative of the 
electron-ion interaction potential appearing in the second 
term of Eq. (86) are given by 



k + G 



d 2 V in 



<9<(q = O)0u?(q= 0) 



k + G' 



(A16) 



-*(G-G')-T S 



-{G a -G' a ){G p -G' p )e 
x (v s (G - G') + <i (k + G, k + G')) • 



$al3h{x) + fl(x)x a X0 



(B2) 



=T S -T,-R 



where the sum over G-space excludes q + G = 0, the 
sums R-space exclude r s — r t — R — 0, and the functions 
fi and f 2 arc defined as follows: 

, 3erfc(^r) + 2 v /^r(3 + 2r,r 2 )e-^ 2 , 
fi(r) = -^-^ — (B3) 



h{r) 



-erfc(^r) - 2yf re^ 2 



(B4) 



APPENDIX B: FORCE CONSTANTS, IONIC TERM 



The total energy of a crystal contains a divergent ion- 
ion energy term, Eq. (3), that combines with the diver- 
gent G = terms in the electron-ion energy (the second 
term in Eq. (18)) to yield the Ewald term Ee w - 



4wN c e 2 ^ e-° / 4 " 
Ee ™ ~ q 2 2^ G 2 

G^O 



5> 



(Bl) 



s,t R 



Z s Z t erfc(v^|T a -T t -R|) 



|r s - r t - R| 



fir; 



E^ 



where erfc(x) = 1 — erf(a;) (erf is the error function); 
the sum over R-space excludes r s — r t — R = 0; A^ c is 
the number of unit cells in the crystal; Z s indicates the 
bare ionic charges for the s-th atom (pseudo-charges in a 
PP-PW framework); 77 is an arbitrary parameter, whose 
value ensures good convergence of both sums over G- and 
R-space. 

The second derivative of the Ewald term yields the 
ionic contribution to the force constants: 

. _ , 47re 2 e -(q+G) 2 /4», 

lon c?f(<i) = E ~(^TGjr ZsZt 

xe^+ G >( T °- T <\q a + G a )(qp + G ) 



2ire 2 



E 



-G 2 /4 V 



fi £-? t G 2 

G^O 

Z s J2Zie lG{T ^G a G p + c.c. 
1 

R 



x 5 af} f 2 {x) + fi{x)x a xp 
e2 ^t EE ZsZl 



x=T 3 -T t -R 



R 
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